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ABSTRACT 


\ 

A  recurring  problem  in  computer  vision  and  related 
fields  is  that  of  generating  computer  models  of  physical  ob¬ 
jects.  This  thesis  presents  a  method  for  constructing  such 
models  in  the  form  of  three-dimensional  surface  or  volume 
descriptions.  The  surface  models  are  composed  of  curved, 
topologically  rectangular,  parametric  patches.  The  data 
required  to  define  these  patches  are  obtained  from  multiple 
photographic  images  of  the  object  illuminated  by  a  rectangu¬ 
lar  pattern  of  lines.  The  projection  of  the  pattern  on  the 
surface  of  the  object  traces  curves  which  define  the 
boundaries  of  the  patches.  The  3D  description  of  the 
patches  is  reconstructed  by  photogrammetric  techniques  from 
two  or  more  images  of  the  projected  pattern.  A  calibration 
stand,  in  which  the  object  is  placed,  permits  determination 
of  the  camera  geometry  directly  from  image  data.—.,  c  * 

This  method  generates  3D  surface  descriptions  of  only 
those  parts  of  the  object  that  are  illuminated  by  the  pro¬ 
jected  pattern,  and  also  are  visible  in  at  least  two  images. 
A  complete  model  of  the  object  is  obtained  by  repeating  this 
reconstruction  process  for  various  arbitrary  orientations  of 
the  object  until  descriptions  covering  the  entire  surface 
have  been  obtained.  Since  each  description  is  defined  in 
its  own  3D  object  space  and  2D  parameter  space,  a  sur¬ 
face-matching  procedure  is  developed  to  register  spatially 


the  surface  descriptions  in  a  common  object  space.  This 
procedure  searches  for  a  3D  rigid  transformation  of  the  sur¬ 
face  descriptions  which  minimizes  their  shape  difference. 
Once  the  surface  descriptions  are  in  the  same  object  space, 
they  are  also  merged  into  a  common  parameter  space.  This 
match-and-merge  process  is  iteratively  repeated  for  pairs  of 
surface  descriptions  until  a  complete  model  of  the  object  is 


assembled 
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CHAPTER  1 


INTRODUCTION 

1.1  Statement  of  the  Problem 

A  recurring  problem  in  computer  vision  and  related 
fields  has  been  automatic  generation  of  computer  models  of 
the  surfaces  of  arbitrarily-shaped,  physical,  three- 
dimensional  objects.  Generation  of  such  models  inherently 
requires  the  acquisition  and  analysis  of  3D  surface  data. 
In  this  context,  acquisition  refers  to  the  ability  to  enter 
automatically  numerical  data  describing  3D  surfaces  into  a 
computer  data  base,  and  analysis  refers  to  the  ability  to 
organize  this  data  base  so  that  it  contains  concise  and 
efficient  models  of  the  complete  surfaces  of  3D  objects. 
The  models  in  such  a  data  base  are  then  available  to  various 
application  tasks.  The  objective  of  the  work  described  here 
was  to  develop  computer  techniques  for  the  construction  and 
processing  of  such  3D  surface  models. 

A  model  of  a  3D  object,  such  as  of  a  machine  part,  can 
be  used  in  a  variety  of  computer  tasks.  In  computer-aided 
design  (CAD),  a  clay  model  of  a  new  part  or  an  existing  part 
is  entered  into  a  computer  data  base  for  engineering  studies 
such  as  finite-element  analysis  (FEA) .  Further  modifica¬ 
tions  to  the  shape  of  the  part  are  made  only  to  the  computer 
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model.  In  ssmputer--a_ided  -manufacturing  (CAM),  the  computer 


model  is  used  to  manufacture  the  part  by  a 
numerically-controlled  machine.  In  robotics,  a  mechanical 
manipulator  with  visual  or  tactile  capabilities  recognizes, 
inspects  and  assembles  the  part  using  the  model.  In 
computer  graphics,  the  model  is  used  for  animated  display  of 
3D  simulation  of  these  tasks. 

In  general,  the  surfaces  being  considered  here  may  vary 
from  those  of  man-made  objects  such  as  industrial  parts  or 
broken  pottery  pieces  to  those  of  natural  objects  such  as 
earth  terrain,  and  the  techniques  developed  here  are 
applicable  to  fields  such  as  automation,  archeological 
restorations,  terrain  recognition,  and  vision  for  intelli¬ 
gent  robots.  This  work  also  attempts  to  bridge  the  gap 
between  the  approaches  to  computer  modeling  of  3D  objects  in 
CAD/CAM  (object  synthesis),  and  those  in  computer  vision 
(object  analysis). 

1.2  Outline  of  the  Approach 

The  process  of  constructing  computer  models  of  surfaces 
of  three-dimensional  objects,  developed  in  this  work,  is 
divided  into  three  major  parts: 

(1)  generation  of  a  numerical  description  of  the  shape 
of  a  surface; 

(2)  hierarchical  structuring  of  this  surface  descrip- 
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tion  to  allow  efficient  processing;  and 
(3)  matching  and  merging  of  partial  surface  descrip¬ 
tions  into  a  complete  surface  model  of  an  object. 

The  surface  representation,  selected  here,  uses  the  bi¬ 
cubic  parametric  patch  as  the  basic  surface  element.  Adja¬ 
cent  patches  with  positional  and  derivative  continuities 
across  their  boundaries  form  sheets  of  composite  patches. 
Given  enough  of  these  composite  patches,  an  arbitrarily- 
shaped  surface  can  be  represented  to  an  arbitrary  precision. 
If  necessary,  however,  this  representation  can  be  expanded 
to  higher  order  bivariate  polynomials,  or  simplified  to 
bilinear  or  planar  patches. 

The  3D  surface  data,  required  to  define  composite  bicu¬ 
bic  patches,  are  obtained  by  photogrammetric  methods  from 
multiple  digital  images  of  surfaces  illuminated  by  a  two- 
dimensional  parameter  space.  The  parameter  space  contains  a 
pattern  of  orthogonal  lines,  defined  on  a  unit-square  grid. 
The  projected  lines  of  the  parameter  space  trace  3D  parame¬ 
tric  curves  on  the  measured  surface.  These  curves  are  the 
boundary  curves  of  the  surface  patches  being  computed.  The 
intersection  points  of  these  curves  are  the  control  points 
of  the  surface  patches.  By  identifying  the  two-dimensional 
projections  of  these  points  and  curves  in  two  or  more  im¬ 
ages,  they  can  be  reconstructed  into  the  original  3D  object 
space.  The  reconstructed  3D  control  points  or  boundary 
curves  are  then  converted  into  linear  lists  of  3D  bicubic 
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parametric  patches. 

Having  generated  a  3D  surface  description  made  of 
composite  patches,  a  hierarchical  control  structure  that 
allows  an  efficient  representation  and  processing  of  these 
patches  is  developed.  Since  the  patches  are  parametrized  in 
a  two-dimensional  parameter  space  where  the  domain  of  each 
patch  is  a  unit  square,  a  hierarchical  structure  in  the  form 
of  a  surface  quadtree  is  utilized.  A  quadtree  of  composite 
patches  is  built  by  recursively  merging  the  four  patches 
that  share  a  common  control  point  into  a  new  patch  at  the 
next  higher  level  of  the  quadtree.  Each  node  in  the  quad¬ 
tree  consists  of  a  description  of  its  patch,  a  bounding 
volume  which  encloses  the  patch,  and  the  average  normal 
vector  of  the  normal  vector  field  of  the  patches  below  the 
current  quadtree  node.  The  bottom  level  of  the  quadtree 
contains  the  original  patches  that  were  generated  from  the 
parametrized  surfaces.  The  higher  levels  contain  coarser 
and  coarser  approximations  to  the  parametrized  surfaces. 
This  quadtree  format  of  bivariate  surface  elements  allows, 
in  general,  logarithmic  rather  than  linear  searching  and 
sorting  time  of  the  individual  elements. 

The  modeling  method  described  up  to  now  generates  only 
a  description  of  a  surface  segment  that  is  parametrized  by  a 
projected  pattern  of  lines,  and  also  visible  in  at  least  two 
images  of  the  surface.  Several  such  segments  must  be  ob¬ 
tained  to  cover  the  complete  surface  of  a  typical  3D  solid 


object.  They  are  obtained  by  changing  the  object's  orienta¬ 
tion  and  repeating  the  surface  acquisition  process.  Each 
surface  segment  is,  therefore,  defined  in  its  own  3D  object 
space.  The  surface  segments  are  obtained,  however,  in  such 
a  way  that  they  partially  overlap  to  allow  their  eventual 
alignment  by  matching  the  common  surface  sections. 

As  the  last  and,  actually,  the  main  step  in  this  mod¬ 
eling  process,  we  match  the  individual  surface  segments  of 
an  object  in  a  common  3D  space,  and  then  also  merge  them 
into  a  common  2D  parametric  space.  A  search  procedure  was 
developed  to  register  spatially  two  partially  overlapping 
surface  segments  at  a  time.  It  employes  a  heuristic-search 
algorithm  with  an  evaluation  function  to  compute  a  rigid  3D 
transformation,  which  minimizes  a  set  of  shape-feature 
distances  between  the  two  surfaces.  The  shape-feature 
distances  can  be  evaluated  at  (1)  patch  control  points,  (2) 
points  with  maximum  surface  curvature,  or  (3)  points 
identified  and  matched  in  images  of  the  surfaces.  These 
feature  points  are  selected  at  each  level  of  a  surface  quad¬ 
tree.  The  search  algorithm  computes  distances  from  the 
feature  points  on  one  surface  to  the  other  surface  by 
tracing  a  ray  normal  to  the  first  surface  at  a  feature  point 
and  intersecting  it  with  the  other  surface.  From  the 
distance,  angle  of  intersection,  and  the  curvature  of  the 
second  surface  at  the  point  of  intersection,  the  algorithm 
evaluates  the  surface  match  of  the  current  node  in  the 


search  tree  and  generates  new  3D  transformations  for  the 
successor  nodes  in  the  search  process.  The  3D  transforma¬ 
tions  of  the  initial  nodes  of  the  search  tree  are  generated 
by  aligning  the  normal  vectors  in  the  top  nodes  of  the  two 
quadtrees  being  matched. 

After  the  two  surface  segments  have  been  transformed 
into  the  same  object  space,  they  must  be  merged  into  a 
common  parameter  space.  Several  algorithms,  which  merge 
overlapping  surface  segments  using  either  transformation  of 
parameters  or  projection  of  parameters,  have  been  developed. 
This  match-and-merge  process  is  iteratively  repeated  for  all 
the  3D  surface  segments  until  a  complete  model  of  the  ob¬ 
ject,  in  a  single  3D  object  space  and  a  single  2D  parameter 
space,  is  generated.  A  surface  model  of  a  solid  object  can 
then  be  converted  into  a  volume  model. 

The  surface-matching  procedure  is  also  generalized  to 
perform  these  two  tasks: 

(1)  surface  and  object  recognition  -  a  partial  3D  sur¬ 
face  description  or  a  complete  object  model  is 
matched  against  a  set  of  complete  object  models  to 
determine  whether  the  surface  may  possibly  be  a 
part  of  one  of  the  objects;  and 

(2)  surface  and  object  segmentation  into  simple  surface 
and  volume  shapes  -  a  3D  surface  description  or  a 
3D  surface  model  of  a  complete  object  is  matched 
with  surface  or  volume  shape  primitives;  once  a 


match  is  obtained,  the  matched  region  is  segmented 
and  the  process  may  again  be  iteratively  repeated 
in  order  to  decompose  a  complicated  object  into  a 
structure  of  simple  shapes. 

The  following  work,  although  not  directly  a  part  of  the 
modeling  process,  has  also  been  performed,  and  is  briefly 
presented  in  this  report: 

(1)  generation  of  synthetic  raster  images  to  allow  a 
realistic  display  of  the  modeled  surfaces; 

(2)  expansion  of  the  traditional  pin-hole  camera  model 
for  display  of  synthetic  raster  images  to  a  model 
of  camera's  optical  system  which  includes  the 
notion  of  focusing,  depth  of  field,  and  motion 
blur  caused  by  a  'lens,  an  aperture  opening,  and  a 
finite  exposure  time,  respectively;  and 

(3)  development  of  the  software  systems  which 
implement  the  modeling  process  as  well  as  the  syn¬ 
thetic  image  generation  and  the  optical  camera 
model. 

Scene  analysis,  image  processing  and  computer  graphics 
are  mostly  applied  sciences.  They  show  that  something  is 
possible  by  doing  it  rather  than  by  formally  proving  that  it 
can  be  done.  The  work  described  here  is  presented  in  this 
spirit  by  ommiting  formal  assumptions,  proofs,  lemmas  or 
anything  else  that  could  just  make  it  look  better  than  it 
actually  is.  It  contains  descriptions  of  techniques,  algo- 


rithms  and  data  structures  and  reports  on  how  well  they  do 
or  do  not  work. 

This  report  has  the  following  structure:  the  three 

steps  of  the  modeling  process  are  presented  in  the  next 
three  chapters;  each  of  these  chapters  is  illustrated  with 
several  examples;  then  a  chapter  is  devoted  to  the  hardware 
employed  and  the  software  developed  for  the  modeling  proc¬ 
ess.  Finally,  a  chapter  with  comprehensive  results  of  the 
modeling  process  follows. 

1.3  Literature  Survey 

There  is  an  extensive  amount  of  recent  literature  that 
is  pertinent  to  this  project  in  the  areas  of  image  proc¬ 
essing,  scene  analysis,  robotics,  artificial  intelligence 
and  computer  graphics.  The  following  brief  survey  is 
divided  into  three  sections  which  correspond  to  the  three 
major  steps  (surface  acquisition,  representation,  and  match¬ 
ing)  in  building  object  models  as  described  in  this  work. 

i*3.i  Sucfece  ftcqtfjsjtiQn 

There  are  a  number  of  techniques  for  automatic,  non¬ 


contact  surface  digitization  based  on  (1)  light  reflection 
from  controlled  light  sources,  (2)  time-of-flight  measure¬ 
ments,  and  (3)  stereometric  correlations.  Illumination  of  a 
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3D  object  by  a  line  pattern  was  first  reported  in  the  image 
processing  literature  by  Will  and  Pennington  [73,74],  who 
used  it  to  locate  planar  faces  of  polyhedral  objects.  A 
normal  vector  to  a  planar  face  'vas  computed  from  a  2D 
Fourier  transform  of  a  single  image  illuminated  by  a  grating 
of  parallel  lines.  Idesawa  et  al  [38,39]  computed  3D  sur¬ 
face  shape  from  Moire  fringe  patterns.  These  patterns  are 
formed  when  an  object  illuminated  by  a  grating  is  observed 
through  another  grating.  Scene  analysis  using  3D  data  ob¬ 
tained  with  a  range  finder  can  be  found  in  several  papers  by 
Agin  [1,2],  Oshima  and  Shirai  [51],  and  Tsukiama  [71]. 
Typically,  a  range  finder  computes  3D  depth  information  from 
a  single  line  of  light  projected  on  a  scene.  Horn  [35]  ob¬ 
tained  surface  geometry  from  single  images  by  analyzing  il¬ 
lumination,  reflectance  and  surface  shading  of  a  scene. 
Posdamer  and  Altschuler  [53]  proposed  a  laser  shutter/space 
encoding  system  capable  of  digitizing  several  thousand  3D 
surface  points  in  real  time. 

Reconstruction  of  objects  and  surfaces  from  multiple 
views  has  been  performed  by  several  authors.  Rabinowitz 
[61]  reconstructed  vertices  and  edges  of  polyhedral  objects, 
Shapira  168]  matched  regions  of  objects  bounded  by  quadric 
surfaces,  Baker  [61  constructed  models  of  objects  from  mul¬ 
tiple  2D  silhouettes.  England  [25]  interactively  fitted  3D 
surface  description  on  two  stereo  views  of  an  object.  Burr 
[13]  matched  two  stereo  images.  Fuchs  et  al  [30]  recon- 
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structed  objects  from  series  of  planar  contours.  A  survey 
by  Bajcsy  [5]  lists  other  current  scene-analysis  methods  for 
acquisition  of  3D  data  using  monocular  and  binocular  depth 


cues. 

1.3.2  Surface  Representation 

Parametric  curved-surface  representation  was  developed 
by  Coons  [20].  The  bicubic  parametric  patch  has  been 
extensively  used  in  computer-aided  design  [7,16,26]  and  in 
computer  graphics  [15,40].  Other  surface  representations 
using  planar  and  quadric  surfaces  were  studied  by  Baumgard 
[10]  and  Levin  [41],  respectively.  Barr  [8]  extended  the 
quadric  representations  into  a  family  of  more  powerful 
shapes.  Volumetric  representation  of  objects  was  used  by 
Agin  and  Binford  [2]  (generalized  cylinders),  O'Rourke  and 
Badler  [50]  (spheres),  and  Meagher  [44]  (octal  trees  of 
cubes).  Hierarchical  surface  representation  using  bounding 
volumes  was  proposed  by  Clark  [18]  to  improve  hidden-surface 
algorithms  for  complex  3D  scenes,  and  implemented  by  Whitted 
[72]  and  Rubin  and  Whitted  [64].  Newell  [47],  Grossman 
[33],  and  Marshall  et  al  [45]  utilize  procedural  representa¬ 
tions  of  3D  objects.  Synthetic  object  modeling  for 
programmed  assembly  by  mechanical  manipulators  can  be  found 
in  works  by  Grossmann  and  Taylor  [32],  and  Lozano-Perez  and 
Winston  [42],  A  survey  of  existing  systems  for  geometric 
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modeling  in  computer-aided  design  is  presented  by  Baer  et  al 
I4J.  Badler  and  Bajcsy  [ 3 ]  survey  methods  of  representing 
3D  objects  for  computer  graphics  and  computer  vision. 

1.3.3  Surface  Matching 

Most  of  the  work  in  surface  matching  and  shape  analysis 
has  been  primarily  done  by  comparison  of  2D  shapes  or  3D 
shapes  projected  into  2D  images.  For  example,  Davis  £ 2 1 J 
used  relaxation  labeling  of  local  features  for  2D  shape 
matching.  Burr  £14]  proposed  a  dynamic  elastic  model  for 
image  as  well  as  line-drawing  matching. 

Relevant  research  in  3D  matching  was  done  by  Baker  £6] 
who  compared  3D  shapes  made  of  piece-wise  surface 
primitives.  Barrow  et  al  £9]  presented  a  technique  of 
matching  a  3D  description  of  a  scene  and  its  2D  images. 
Horn  and  Bachman  £36]  aligned  a  synthetically  shaded  image 
of  a  3D  terrain  model  with  an  actual  image  of  the  terrain. 
Burr  £131  matched  reconstructed  3D  edges  of  objects  with 
stored  wire  frame  models  utilizing  3D  features  and  geome¬ 
trical  constraints.  Shapiro  et  al  £691  matched  structured 
descriptions  of  objects. 

1.4  Summary  of  Definitions 

The  modeling  process  described  here  utilizes  three  co¬ 
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ordinate  systems  and  four  mappings  between  these  systems. 
We  need  to  define  these  systems  and  mappings  in  detail 
before  proceeding  further: 


0(x/y,z,w)|?  is  a  homogeneous  3D  object  coordinate  sys¬ 
tem  which  defines  object  space  k.  All  3D  representations  of 
surface  elements  are  specified  in  an  object  coordinate  sys¬ 
tem.  An  object  being  modeled  is  positioned  in  object  coor¬ 
dinate  spaces  Ji  =  1/  2,  ....,  J&.  Notation  0  (x,y,  z,w)  ^  ^ 
means  that  object  space  j£  is  parametrized  by  parameter  space 

j- 

P(urvfw)j  is  a  homogeneous  2D  parameter  coordinate  sys¬ 
tem  which  defines  parameter  space  j.  A  set  of  orthogonal 
lines  on  a  unit-square  grid  parametrizes  a  3D  surface  in  an 
object  space.  An  object  being  modeled  is  parametrized  by 
parameter  coordinate  systems  j  =  1,2,  ....,  J[.  There  is 
always  at  least  one  parameter  coordinate  system  j  for  every 
object  coordinate  system  Jc.  Notation  PtUrVrw).^^  means  that 
parameter  space  j  parametrizes  object  space  k;. 

Q(rfsfw)i  is  a  homogeneous  2P  imgge _ SggLdijiate _ syst-SlP 

which  defines  image  space  i.  An  object  being  modeled  is 
observed  in  image  coordinate  systems  i  =  1,  2,  ....,  1. 
There  are  always  at  least  two  image  coordinate  systems  i  and 
i’  for  every  object  coordinate  system  Jc  and  parameter  coor¬ 
dinate  system  j.  Notation  Q (r, s, w) ^ f j f k  means  that  object 
space  Js  is  parametrized  by  parameter  space  j  and  observed  in 
image  space  i. 


T ( P j — >0^)  is  a  mapping  from  a  2D  parameter  space  Pj  to 
a  3D  object  space  Ok.  This  mapping  takes  places  (via  a  pro¬ 
jector)  when  2D  lines  of  a  parameter  space  are  projected  on 
a  3D  object  space. 

TfO^— )  is  a  mapping  from  a  3D  object  space  0k  to  a 
2D  image  space  Q^.  This  mapping  takes  place  (via  a  camera) 
when  a  2D  projection  of  a  parametrized  surface  is  recorded 
in  a  2D  image  coordinate  system. 

TfO^-^O^i)  is  a  mapping  from  object  space  0^  to  object 
space  0ki.  This  mapping  is  computed  by  the  matching  proce¬ 
dure  to  transform  two  3D  surface  descriptions  into  a  common 
object  space. 

T{Pj-»Pji)  is  a  mapping  from  parameter  space  Pj  to  pa¬ 
rameter  space  Pj i .  This  mapping  is  computed  by  the  merge 
procedure  to  transform  a  3D  surface,  parametrized  in  two 
different  parameter  spaces,  into  a  common  parameter  space. 

Additionally,  the  following  terms  are  used  throughout 
this  report: 

A  parametrizing  grid  is  a  set  of  orthogonal  2D  isopa¬ 
rametric  lines  used  to  parametrize  a  3D  surface. 

An  isoparametric  line  network  is  a  2D  image  projection 
of  a  parametrizing  grid  projected  on  a  3D  surface. 

A  surface  element  is  a  3D  surface  representation  by  a 
planar  face,  a  sphere,  a  quadric  surface,  or  a  bicubic 
patch. 

A  bounding  surface  is  a  surface  element  used  solely  as 


a  bound  of  other  surface  element  or  elements. 


A  bounding  volume  is  a  set  of  one  or  more  surface  ele¬ 
ments  which  completely  enclose  a  3D  space.  One  or  more  sur¬ 
face  elements  or  other  bounding  volumes  are  completely 
circumscribed  by  a  bounding  volume.  Typical  bounding 
volumes  are  a  sphere,  an  ellipsoid  or  a  parallelepiped. 

A  surface  quadtree  is  a  hierarchical  surface  structure 
with  surface  representation  parametrized  by  two  orthogonal 
parameters  which  are  defined  in  a  2D  plane.  A  node  of  the 
tree  has  up  to  four  successors  which  contain  descriptions  of 
the  current  surface  element  subdivided  into  four  more 
detailed  elements. 

A  sheet  of  composite  patches  is  a  set  of  contiguous 
patches  with  adjacent  patches  having  at  least  first- 
derivative  continuity  (C*)  across  their  boundaries.  All 
patches  in  a  sheet  are  parametrized  by  the  same  parameter 
coordinate  system.  A  sheet  of  patches  is  stored  in  a  single 
surface  quadtree.  Patches  in  different  sheets  of  an  object 
have  at  most  positional  continuity  (C®) . 

A  surface  segment  is  the  part  of  a  surface  that  is 
described  in  one  object  coordinate  system  and  one  parameter 
coordinate  system. 

A  surface  model  of  a  solid  object  is  a  union  of  sheets 
of  composite  patches  which  define  the  surface  of  the  object. 

A  volume  model  of  a  solid  object  is  a  union  of  rectan¬ 
gular  parallelepipeds  which  define  the  volume  of  the  object. 
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CHAPTER  2 


GENERATION  OF  3D  SURFACE  DESCRIPTIONS 

This  chapter  describes  a  method  for  making  automatic, 
non-contact  measurements  of  the  surfaces  of  physical,  rigid, 
arbitrarily-shaped  3D  objects  [54,28,56].  Surface  informa¬ 
tion  is  obtained  by  photogrammetric  and  image-processing 
techniques  from  multiple  images  of  the  measured  surfaces  il¬ 
luminated  by  a  controlled  light  source.  The  method,  as  used 
here,  generates  surface  data  that  are  used  in  modeling  the 
surfaces  with  composite  bicubic  patches.  However,  it  can 
also  be  adapted  to  surface  modeling  with  surface  elements 
ranging  from  planar  triangles  to  high-order  polynomials. 

Surface  representation  by  the  bicubic  parametric  patch, 
as  defined  in  detail  in  Appendix  A,  has  been  chosen  because 
of  its  widespread  applications  in  the  fields  where  computer¬ 
generated  models  of  3D  objects  are  extensively  used,  such  as 
CAD/CAM  and  computer  graphics.  The  format  of  the  patch 
employed  here  can  be  defined  either  by  a  set  of  3D  control 
points  or  by  a  set  of  3D  boundary  curves.  The  photo¬ 
grammetric  reconstruction  method  allows  us  to  compute  either 
the  control  points  or  the  boundary  curves  from  multiple  im¬ 
ages  of  the  surface-patch  area. 

The  first  three  sections  of  this  chapter  contain  the 
background  mathematics  necessary  for  the  reconstruction 


process.  The  fourth  section  then  describes  the  actual 
computer  methods  implementing  this  process. 


To  model  an  existing  3D  surface  with  parametric 
patches#  the  surface  is  parametrized  with  orthogonal  isopa¬ 
rametric  lines.  T(Pj— >0^}  maps  a  2D  parameter  homogeneous 
coordinate  system  P(u#v#w)j  into  a  3D  object  homogeneous  co¬ 
ordinate  system  Otx^yfZ^w)^  [Figure  2.1]: 

(2.1) 


The  parametric  information  mapped  on  the  surfaces  in 
the  object  space  can  be  described  by  a  set  of  orthogonal 
lines  and  their  intersection  points.  There  are  (15+1)  x 
(JS+1)  isoparametric  lines;  those  with  u  as  the  parameter 
and  y  constant  are  labeled: 

1q(u)#  . . .  r  1  n  (u ) '  *  •  *  t  ljj(u)# 

and  those  with  y  as  the  parameter  and  jj  constant  are 
labeled: 

T0(v),  ...,  Tm(v),  ...,  Tm(v). 

The  intersection  point  of  two  orthogonal  lines  Tn(u)  and 
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Tm(v)  is  a  patch  control  point  T^^tu.v.w)  [Figure  A. 1(b)]. 
This  point  is  mapped  on  surface  point  n  (x,  y,  z, w)  k  and 
then  projected  into  image  point  h^ n (r, s,w) ^ .  It  becomes  a 
control  point  of  patch  gmrn(u,v).  The  line  segment  from 
point  ^m,n(u'v'w)  to  point  ^>,,*+1 (u' v'w>  islm^n(v).  It  is 
mapped  into  a  3D  patch  boundary  curve  cmfn(v)k  =  [  x(v)  y(v) 
z(v)  w(v)  ]k  and  then  projected  into  a  2D  image  curve 
cm,n^i  s  ^  c(v)  s(v)  w(v)  1^.  The  other  line  segment  from 
point  to  point  Vn,  n  v,  w)  is  similarly  label¬ 

ed.  The  mappings  of  these  two  line  segments  become  boundary 
curves  of  patch  gm^n(u,v).  This  defines  the  surface  infor¬ 
mation  required  to  construct  an  array  (a  sheet)  of  M  x  N  pa¬ 
rametric  patches. 

2.2  3D  to  2D  Image  Mapping 


Transformation  T(0k— maps  a  3D  object  homogeneous 
coordinate  system  0(x,y,zfw)k  to  a  2D  image  homogeneous  co¬ 
ordinate  system  QUfSrW^  [Figure  2.21: 


(2.2) 
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The  transformation  T{Ok-*^Q^}  consists  of  translations 


and  rotations  of  the  camera  in  0 (x, y, z , w) k,  perspective  pro¬ 
jection,  and  conversion  into  the  digitized  image  coordinate 


system  Q(r,s,w)^.  Its  detailed  composition  has  been 

previously  discussed  in  [54]. 

2.3  2D  to  3D  Surface  Reconstruction 

To  map  surface  information  from  an  image  plane 
Q(r , s,w) ^  back  into  the  3D  space  0(x,y,z,w)k  we  have  to  find 
an  inverse  mapping  to  that  given  in  equation  (2.2).  For  the 
2D  projection  Ti(r,s,w)^  in  Q(r,s,w)^  of  a  3D  point 

h(x,y,z,w)k  in  0(x,y,z,w)k  and  mapping  from  the  3D 

object  coordinate  system  to  the  2D  image  coordinate  system, 
we  can  rewrite  equation  (2.2)  into  a  system  of  two  linear 
equations  with  three  unknowns: 

(2.3) 

rit32"t12  rit33"fc13  rit34"t14 

sifc32_t22  sit33~t23  sit34“t24_ 


The  solution  to  this  system  of  equations  produces  only 
the  equation  of  the  3D  line  H(r,s,w)^  h (x, y, z, w) k.  To 
compute  the  coordinates  of  the  point  "E(x,y,z,w)k  we  need  an 
additional  projection  E(r,s,w)^i  in  Q(r,s,w)^i  and  mapping 
T(Ok-»Q^i}  as  shown  in  Figure  2.2.  With  two  projections  we 
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obtain  four  equations  with  three  unknowns  and  the 
coordinates  of  point  E (x, y, z, 1) ^  can  be  solved.  The 
least-square  error  solution  of  the  system  in  equation  (2.3) 
for  1  projections  was  previously  described  in  [54],  The 
reconstruction  of  surface  points,  presented  in  this  section, 
can  be  generalized  to  reconstruction  of  surface  curves. 
Appendix  B  contains  an  analogous  method  for  the 
reconstruction  of  a  3D  parametric  cubic  curve  from  its 
multiple  projections  in  2D  images. 

The  reconstruction  error  of  a  point  E(x,y,z,l)|.  is 
evaluated  from  equation  (2.2)  by  projecting  the 
reconstructed  point  back  into  each  image  i  from  which  it  was 
obtained,  and  computing  the  2D  image  distance  between  this 
projection  and  the  originally  measured  point  E (rif si , 1) A . 
The  total  error  is  then  the  sum  of  errors  in  the  1  images: 


,  tiix+ti2y+ti3z+ti4  t21x+t22y+t23z+t24,? 

Jss2-Yri - >  +  {si - > 

i=l  t31x+t32y+t33z+t34  fc31x+t32y+ fc33z+ fc34 


In  the  above  equation  (2.4)  the  coefficients  t  belong 
to  the  camera  transformation  matrix  T(0k— »Q^}.  Note  also 
that  the  error  d  is  given  in  pixel  units  of  the  image 
coordinate  systems. 

The  twelve  coefficients  of  the  camera  transformation 
matrix  TCO^— in  equation  (2.3)  are  computed  for  each 


age  i.  by  a  camera  calibration  method  [703.  Equation  (2.2) 
may  be  rewritten  as  two  simultaneous  linear  equations  by 
eliminating 

(2.5) 

x  y  z  1  0  0  0  0  -xr  -yr  -zr  -r~|  I"  tn 

OOOOxyzl  -xs  -ys  -zs  -s  ti? 

•  •  • 

•  •  • 

fc33 
_  fc34 

where  x,  .y,  z,  are  the  coordinates  of  a  point  Ti(x,y,z,l)k  in 
0(x,y,z,w)k  and  X/  £  ate  the  coordinates  of  its  image 
h(r,s,l)^  in  Q(r,s,w) There  are  twelve  unknown  i's  in  the 
two  equations  of  (2.5).  If  we  know  six  points  in 

0(x,y,z,w)k  and  their  corresponding  images  in  Q(r,s,w)i,  we 
can  create  twelve  simultaneous  equations  (two  for  each 
point)  and  solve  for  the  unknown  coefficients.  In  order  to 
solve  such  a  system  of  linear  equations,  it  has  to  be  made 
non-homogeneous  be  letting,  for  example,  t34  =  1.  This  is 
allowed  since  the  homogeneous  coordinate  systems  that  we  are 
using  have  an  arbitrary  scale  that  is  set  to  1. 

In  general,  for  six  or  more  calibration  points  the 
least-squares  error  method  is  again  used  to  solve  a  system 
of  eleven  non-homogeneous  equations  (with  t34  =  1)  [543. 

Once  the  coefficients  of  T(Ok-»Q^}  are  computed,  the  error 
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of  each  calibration  point  may  be  checked  by  evaluating 
equation  (2.2)  for  h(x,y,z,l)k  and  comparing  the  computed 
point  h(r,s,l)^  with  the  measured  calibration  point 
h(r,s,l)^.  Calibration  points  with  large  errors  may  then  be 
disregarded  and  the  procedure  repeated  for  a  better  estimate 
of  T(Ok-»Qi}. 

2.4  Photoorammetric  Reconstruction 

In  the  approach,  taken  here,  the  surface  to  be 
digitized  is  illuminated  by  the  parametrizing  pattern 
described  in  Section  2.1.  The  pattern  is  focused  on  the  3D 
surface  and  its  shadow  traces  3D  isoparametric  curves  there. 
Using  image  processing  and  photogrammetric  techniques,  a  3D 
description  of  the  pattern  projected  on  the  surface  can  be 
obtained.  Therefore,  the  3D  surface  being  measured  is 
illuminated  by  a  2D  pattern  of  parametric  lines  P(u,v,w)j 
with  a  projector  located  at  j  (Figure  2.11.  The  illuminated 
parts  of  the  surface  are  photographed  from  at  least  two 
locations  i  and  i'  by  a  camera  onto  2D  image  planes 
Q(r,s,w>£  and  Q(r,s,w)^«,  respectively  (Figure  2.2).  The  3D 
object  coordinate  system  0(x,y,z,w)k  is  defined  by  a  set  of 
camera  calibration  marks  placed  in  a  camera  calibration 
stand.  The  measured  surface  is  within  this  stand. 

This  arrangement  limits  us  to  surfaces  that  can  be  pho¬ 
tographed  under  the  conditions  imposed  by  this  method.  That 
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is,  the  surfaces  are  photographed  in  a  camera  calibration 
stand;  they  can  not  be  transparent,  have  high  specular  re¬ 
flection,  or  be  totally  black.  This  method  also  requires  a 
large  amount  of  computations,  and  image  digitization  at  a 
high  resolution  to  obtain  precise  surface  measurements.  It 
is  mainly  intended  for  surface  measurements  of  objects  w^ose 
models  will  later  be  used  in  various  applications,  rather 
than  for  real-time  surface  acquisition  and  processing. 

The  following  is  an  outline  (specified  in  a  home-grown 
variation  of  the  "Algol-like  notation")  of  the  recon¬ 
struction  algorithms,  once  all  digitized  images  for  a 
single  parameter  coordinate  system  j  and  a  single  object  co¬ 
ordinate  system  k.  have  been  obtained: 

(2.6) 


procedure  extract  (I j f ^ ) ; 

begin 

l£LL  i  :=  1  Step  1  until  I-i  do. 
begin 

extract  coordinates  of  camera  calibration  marks; 

match  them  with  their  coordinates  in  0(x,y,z,w)^; 

compute  TtO^-^Q^}; 

extract  the  projected  line  pattern; 

end; 

E.eturn; 

end; 


and  then 

procedure  reconstruct  (1^  j^); 

begin  J' 

match  extracted  patterns  for  images  1  to  1^  j.; 
reconstruct  matched  surface  data  for  images  I  to 
return; 
end; 


(2.7) 


Images  of  surfaces  are  taken  by  a  camera  on  monochrome 
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film  and  digitized.  There  are  two  types  of  data  that  need 
to  be  extracted  from  a  digitized  image:  the  projected  line 
patterns  and  the  illuminated  calibration  marks.  Both  data 
types  can  be  stored  in  binary  images  (each  pixel  is  either  0 
or  1),  where  the  projected  pattern  is  a  network  of  black 
lines  on  a  white  surface  with  black  background  (Figure 
2.3(a)],  and  the  calibration  marks  are  white  marks  on  black 
background  (Figure  2.3(b)]. 

The  projector  illuminating  the  measured  surface  can  be 
considered  a  single  point-light  source.  The  light  reflec¬ 
tion  from  the  surface  is  only  diffuse  reflection;  any 
potential  specular  reflection  can  be  supressed  by  the  use  of 
polarizing  filters.  Therefore  the  light  intensity  in  an  im¬ 
age  can  be  approximately  modeled  by  Lambert's  law: 

intensity  =  k^  (N*H)  (2.8) 

where 

k^  =  diffuse  reflection  coefficient 

N  =  surface  normal  vector 

H  =  light  source  vector 

intensity  =  final  image  intensity 

It  is  apparent  from  equation  (2.8)  that  image  intensity 
and,  therefore,  also  the  contrast  of  the  projected  lines 
decrease  as  the  illuminated  surface  turns  away  from  the  pro¬ 
jector.  Simultaneously,  the  spacing  of  the  lines  projected 
into  the  image  decreases  as  does  their  width.  A  frequency- 


domain  filter,  or  edge  detection  are  used  to  restore 
contrast  to  these  parts  of  an  image. 

Each  image  is  digitized  to  intensity  resolution  of  8  or 
16  bits  per  pixel.  The  digital  image  is  then  filtered  to 
remove  noise  and  to  enhance  contrast  before  being 
thresholded  into  a  binary  image.  Image  noise  is  removed  by 
spatial  convolution  with  a  3x3  or  5x5  pixel  operator. 
Contrast  enhancement  of  the  line  pattern  is  performed  by  (a) 
frequency-domain  filter,  (b)  local  histogram  modifications 
in  small  image  windows,  or  (c)  edge-detecting  heuristic 
search  [43]  which  is  guided  by  the  shape  of  the  projected 
lines  and  the  tendency  of  adjacent  lines  to  follow  similar 
shape.  The  frequency-domain  filter  is  a  band-emphasis  and 
low- pass  filter  [34,60]  designed  to  enhance  the  frequencies 
of  the  projected  grid  lines,  and  to  remove  high-frequency 
noise. 

A  program  computes  the  coordinates  of  all  the  visible 
corners  of  camera  calibration  marks  by  taking  these  steps  on 
a  binary  image: 

(2.9) 

(1)  read  a  window  of  the  binary  image; 

(2)  make  a  fast  check  for  a  possible  presence  of  a 
calibration  mark  in  the  window,  if  not  present  go 
to  (1);  (this  fast  check  is  made  by  parsing  pixels 
in  the  x  and  x  directions  and  looking  for  pixel 
strings  of  types  0raln,  lm0n,  or  lm0nlP; 


(3)  compute  edge  pixels  of  a  calibration  mark  with  an 
8-point  Laplacian  operator; 

(4)  locate  straight  lines  in  edge  pixels  using  Hough 
transform  [23]; 

(5)  compute  parametric  equations  of  lines  fitted  into 
these  edge-lines; 

(6)  compute  intersections  of  the  lines;  these  2D  in¬ 
tersection  points#  h(r,s)^,  are  the  corners  of  a 
camera  calibration  mark#  go  to  (1) . 

Next#  each  extracted  corner  (r#s)  is  matched  with  its 
pre-stored  coordinates  in  0(x,y#z,w)k  and  a  five-tuplet 
(r#s;x,y#z)  of  the  coordinates  is  formed.  The  camera  trans¬ 
formation  matrix  T{Ok—»Q^)  is  then  computed  from  six  or  more 
of  these  five-tuplets  by  solving  the  system  of  equations 
(2.5).  Since  the  calibration  stand  is  normally  photographed 
by  the  camera  only  from  one  octant  of  the  object  coordinate 
space  0k,  a  simple  relational  graph  [ 2 4 # 75 ]  has  been  devel¬ 
oped  to  perform  the  matching  of  the  projected  calibration 
points  with  their  3D  coordinates.  The  graph  contains  2D 
spatial  relationships  of  the  calibration  marks  when  viewed 
by  the  camera  from  the  given  octant. 

The  networks  of  the  projected  parametric  lines  are  ex¬ 
tracted  from  a  binary  image  by  a  program  which  may  execute 
concurrently  on  the  same  image  with  the  program  extracting 
the  camera  calibration  marks.  This  program  creates  a 


network  of  lines  by  performing  these  steps: 

(2.10) 

(1)  apply  a  line- thinning  operator  to  the  image;  this 
operator,  similar  to  those  given  in  [52,763,  uses  a 
transition  table  to  thin  iteratively  lines  in  a 
binary  image  while  preserving  line  connectivity, 
i.e.,  a  line  pixel  has  only  1  or  2  adjacent  pixels 
while  a  node  pixel  has  3  or  4  adjacent  pixels  in  an 
8-pixel  neighborhood; 

(2)  remove  all  pixels  which  do  not  belong  to  lines  or 
nodes,  i.e.,  pixels  with  0  or  more  than  4  adjacent 
pixels  (background,  calibration  marks,  shadows,  and 
surfaces  not  illuminated); 

(3)  convert  the  line  pixels,  still  in  raster  format, 

into  parametricly  defined  straight  line  segments  or 
cubic  curves;  this  is  the  scan-line-to-vector 
conversion  where  all  adjacent  line  pixels, 

connecting  two  node  pixels,  are  collected  into 
clusters  of  ”E(r,s)^  points  and  either  a  2D  straight 
line  segment  or  a  cubic  curve  is  fitted  into  them, 
the  two  end  points  of  the  line  segment  or  the  curve 
are  also  computed. 

The  above  three  steps  are  performed  in  one  pass  over  a 
binary  image  while  no  more  than  10  to  15  scan  lines  have  to 
be  accessible  at  the  same  time.  After  the  image  has  been 
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processed,  node  intersections  are  computed.  For  the  pattern 
of  orthogonal  parametric  lines  used  here,  there  are 
typically  four  line  segments  or  curves  intersecting  at  a 
node.  The  coordinates  of  a  node  are  computed  by  minimizing 
the  sum  of  distances  to  the  lines  or  curves  from  the  inter¬ 
section  point.  A  data  structure  is  created  for  the  line 
network  where  each  node  description  contains  the  h(r,s)  co¬ 
ordinates  of  the  node,  pointers  to  all  the  adjacent  nodes 
and  pointers  to  coefficients  of  curves  or  lines  connecting 
adjacent  nodes.  Pointers  of  each  node  to  its  neigboring 
nodes  are  sorted  in  such  a  way  that  each  isoparametric  line 
forms  a  two-way  linked  list  of  its  segments.  Therefore  the 
topology  of  the  line  network  is  implicitely  embedded  in  the 
data  structure  and  one  may  traverse  the  line  network  in  any 
path.  This  capability  is  used  during  the  reconstruction  of 
the  projected  lines. 

It  should  be  emphasized  that  while  computing  the  2D 
projections  of  the  surface  curves,  there  is  additional  in¬ 
formation,  not  needed  for  the  reconstruction  of  individual 
3D  points,  that  has  to  be  extracted  from  each  image.  When 
constructing  cubic  curves  with  the  method  of  Appendix  B,  the 
slope  of  the  projected  curve  at  each  end  point  must  also  be 
evaluated. 

Finally,  a  matching  algorithm  which  reconstructs  the 
multiple  2D  projections  of  the  isoparametric  line  networks 
in  the  3D  object  coordinate  system  is  presented.  This  algo- 
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rithm  matches  corresponding  projections  of  points  and  carves 
which  are  then  numerically  reconstructed  by  equation  (2.3) 
for  points,  and  the  method  of  Appendix  B  for  curves.  This 
matching  is  primarily  based  on  the  topology  of  the  projected 
line  networks.  The  matching  proceeds  in  two  steps:  (1)  an 
initial  estimate  of  the  match  is  made,  and  (2)  the  correct 
match  is  found  by  a  gradient  search. 

In  the  first  part,  the  algorithm  locates  a  cycle  of 
four  adjacent  nodes  in  one  image  and  attempts  to  find  pro¬ 
jections  of  these  nodes  in  the  other  image.  The  projections 
must  also  form  a  cycle  of  adjacent  nodes.  The  projection  of 
a  node  is  determined  when  the  3D  line  from  the  center  of 
projection  to  the  node  in  the  first  image  is  projected  into 
the  second  image  and  the  node  nearest  to  it  is  found.  This 
matches  node  n ( r , s ) ^  with  node  hm^  n  (r ,  s)  ^  ,  and  similarly 
the  other  three  nodes  [Figure  2.4(a)].  Note  that  the  four 
nodes  in  a  cycle  provide  the  minimum  surface  information 
needed  for  a  single  patch  element.  This  initial  match 
succeeds  if  the  reconstruction  error,  evaluated  by  equation 
(2.4),  of  each  of  the  four  nodes  is  less  than  a  specified 
maximum.  The  acceptable  maximum  error  depends  on  the  reso¬ 
lution  of  the  digitized  images.  If  it  fails  then  another 
cycle  is  located  and  this  step  is  repeated;  if  all  the 
cycles  are  exhausted  without  a  success,  then  the  two  pro¬ 
jections  do  not  correspond  to  each  other  (either  they  con¬ 
tain  different  surfaces,  or  they  are  parametrized  by  differ- 


other  (either  they  contain  different  surfaces,  or  they  are 
parametrized  by  different  parameter  systems).  Although  we 
could  match  only  a  single  node  in  this  initial  step,  a  cycle 
of  four  nodes  substantially  improves  the  match  and, 
therefore,  shortens  the  search  in  the  second  step. 

In  the  second  part,  an  evaluation  function  is  used  to 
search  for  the  correct  match  of  the  two  line  networks. 
First,  the  current  match  of  nodes  Tim^n(r,s)i  and  hm^n(r,s)it 
is  iteratively  propagated  over  the  line  network  in  the 
directions  of  the  parametric  lines  until  all  matchable  nodes 
are  reconstructed.  Then  the  evaluation  function  which 
minimizes  the  sum  of  the  node  reconstruction  errors  and 

maximizes  the  number  of  reconstructed  nodes  is  computed  by: 

L 

e(m,n;m,n)  *  wdSdp  +  Wl/L  (2.11) 

P=1 

where 

dp  =  reconstruction  error  of  node  £ 

wd  =  weight  of  ^dp 

L  =  number  of  reconstructed  nodes 

wL  =  weight  of  1/L 

e  =  value  of  evaluation  function  for  match  of  node 
Kfnl*'s)i  with  node  ^m,  n  ^r'  s)  i  1  • 

The  first  line  network  is  then  moved  by  one  node  in  the 
four  directions  of  the  network  with  respect  to  the  second 
network  and  four  new  matches  of  node  n  (r ,  s)  ^  with  nodes 


i 
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The  evaluation  func- 


s ) ^ i  are  generated  [Figure  2.4(b)]. 
tions  e  (m,  n;m+l, n) ,  e  (m,  n;m-l,  n) ,  e  (m,  n;m,  n+1) ,  and  e(m,n; 
m,n-l)  of  each  match  are  computed,  and  the  correct  match  is 
found  when  the  value  of  the  function  for  a  match  is  smaller 
than  the  value  at  any  of  its  neighbors,  that  is: 

e(m,n;m,n)  <  min£  e  (m,n;m+l,n)  ,e  (m,n;m-l,n) ,  (2.12) 

e  (m,  n;m,  n+1) , e  (m,  n;m,  n-1)  ~] 

If  the  inequality  (2.12)  is  not  satisfied  by  the 
current  match  then  the  match  with  the  smallest  value  of  e 
becomes  the  current  match  and  this  search  is  repeated  until 
it  converges  on  the  solution  and  the  inequality  (2.12) 
becomes  valid.  The  final  match  is  then  used  to  create  a 
data  structure  which  contains  topological  description  of  the 
reconstructed  parts  of  the  line  network  again  defined  as 
two-way  linked  lists  of  the  orthogonal  line  segments  and  the 
reconstructed  3D  points  and  curves.  This  data  structure  can 
be  converted  into  a  linear  list  of  surface  data  needed  to 
compute  individual  bicubic  patches. 

If  the  images  contain  two  or  more  disconnected  line 
networks,  each  network  has  to  be  processed  by  this  algorithm 
separately.  For  more  than  two  projections,  this  process  is 
carried  out  for  the  first  two  projections  as  described  and 
then  each  additional  projection  is  matched,  one  at  a  time, 
with  the  already  reconstructed  network. 


In  many  applications  of  surface  models,  specially  mod¬ 
els  of  closed  (solid)  objects,  it  is  pertinent  to 

distinguish  between  the  "outside"  and  the  "inside"  sides  of 
a  surface.  The  "outside  of  a  surface  or  object  is  defined 
as  the  side  facing  the  cameras  during  the  surface  recon¬ 
struction.  A  surface-normal  vector,  pointing  to  the 
"outside"  direction,  is  attached  to  each  reconstructed  node. 
The  vector  is,  therefore,  forced  to  be  oriented  so  that: 

N  •  Hj.  >  0  (2.13) 

where 

N  ■  surface-normal  vector. 

Hi  =  location  of  center  of  projection  of  image  i. 

The  center  of  projection  of  image  i  is  computed  from  the 

system  of  equations  (2.3)  using  the  camera  transformation 
matrix  TCO^— and  two  linearly  independent  points  in  the 
image  plane.  Each  equation  of  (2.3)  defines  a  3D  plane; 
their  intersection  is  a  3D  line  through  the  center  of  pro¬ 
jection.  The  intersection  of  two  such  lines  is  the  location 
of  the  center  of  projection,  denoted  H^.  If  only  two  images 
are  used  to  generate  the  3D  data,  then  either  center  of  pro¬ 
jection  can  be  used  in  (2.13)  to  orient  the  normal  vector 
since  each  point  must  have  been  visible  in  both  images 
(otherwise  it  could  not  have  been  reconstructed).  If  more 
than  two  images  are  used,  then  each  point  needs  a  pointer  to 

the  camera  transformation  matrix  of  one  of  the  images  from 


which  it  was  reconstructed.  The  location  of  the  center  of 


projection  of  this  image  is  used  in  (2.13). 

The  surface  information,  computed  at  each  reconstructed 
node,  consists  of  the  following  items: 

h(x,y,z)  positional  coordinates  (2.14) 

surface  tangent  in  u 
surface  tangent  in  v 
surface  cross-derivative  (twist) 
surface-normal  vector 

pointer  to  adjacent  node  in  +u  direction 
pointer  to  adjacent  node  in  +v  direction 
pointer  to  adjacent  node  in  -u  direction 
pointer  to  adjacent  node  in  -v  direction 
The  above  information  is  used  to  convert  the  reconstructed 
nodes  and  curves  into  surface  patches  as  outlined  in 
Appendix  A. 


hu  (x,y,  z) 
hv (x, y, z) 
huv(x,y,z) 
N (x,y, z) 
pointer^ 
pointer 2 
pointer3 
pointer4 


2.5  Illustrations 


The  surface  reconstruction  method  is  illustrated  here 
by  the  following  example  of  a  multiply-curved  surface. 
There  are  two  images,  in  Q(r,s,w)1^1^1  and  Q(r,s,w)2f i  i  im¬ 
age  coordinate  systems,  of  the  parametrized  surface  as  shown 
in  Figure  2.5.  The  images  are  displayed  at  resolution  of 
1100  x  1400  binary  pixels.  The  orthogonal  parametric  lines 
extracted  from  the  two  images  are  drawn,  together  with  a 
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model  of  the  camera  calibration  marks  and  the  object  coordi¬ 
nate  system  0<x,y,z,w)2,  in  Figure  2.6.  The  image  recon¬ 
struction  mappings  TtO^-^Q^}  and  T{0^— >Q2}  are  also  printed 
in  this  figure  under  their  corresponding  images.  The  recon¬ 
structed  surface,  named  'SURFACE. 3',  is  shown  in  four  ortho¬ 
graphic  projections  in  Figure  2.7.  The  surface  is 
represented  by  a  network  of  3D  polygons.  There  are  four 
different  projections  of  the  surface  shown  in  Figure  2.8. 
Here,  the  surface  is  represented  by  a  sheet  of  composite  3D 
bicubic  patches  with  C*  continuity.  Shaded  synthetic  images 
of  the  sheet  of  patches  are  shown  in  Figure  2.9. 

2.6  Summary 

A  technique  for  reconstruction  of  3D  surface  informa¬ 
tion  from  multiple  images  was  presented  in  this  chapter.  It 
is  based  on  correlation  of  projected  patterns  of  light  on 
the  surfaces  being  measured.  Triangulation  is  used  to 
compute  3D  surface  data  from  2D  data  obtained  in  two  or  more 
images.  Image  transformations  are  computed  from  a  number  of 
calibration  marks  whose  locations  in  the  3D  space  are  known, 
and  whose  locations  in  an  image  are  measured.  Arbitrary 
camera  positions  and  orientations  are  therefore  permitted, 
however,  limitations  of  stereo  correlation  still  persist: 
reconstruction  error  increases  with  narrow  separation  angle, 
reconstruction  data  decrease  with  wide  separation  angle. 
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This  method  allows  additional  images  to  be  used  to  improve 
accuracy  of  the  measured  surfaces  and  to  increase  the 
measured  surface  area.  However,  all  surface  measurements  in 
this  work  were  obtained  from  pairs  of  images  sequentially 
taken  by  the  same  camera. 
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CHAPTER  3 


HIERARCHICAL  SURFACE  REPRESENTATIONS 

Computer  surface  modeling  refers  to  the  ability  to 
analyze  a  three-dimensional  surface  and  represent  it  by  a 
composite  description  based  on  many  independent  features 
such  as  size,  shape,  structure,  texture,  and  light  reflec¬ 
tive  and  transmissive  characteristics.  The  creation  of  such 
a  computer  model  is  necessary  for  any  computer  vision  or 
scene  analysis  system.  This  chapter  describes  a  surface  and 
volume  representation  system  which  was  developed  to  support 
the  modeling  process  of  this  work. 

The  representation  of  3D  surfaces  in  a  data  base  is  hi¬ 
erarchical?  that  is,  it  consists  of  a  control  structure 
which  in  the  most  general  case  is  a  directed,  acyclic  graph. 
The  nodes  of  the  graph  contain  descriptions  of  the  surface 
shapes  and  their  relations,  and  the  arcs  represent  3D  rigid 
transformations  in  the  object  coordinate  system.  There  are 
three  types  of  nodes  in  such  a  graph:  (1)  surface  elements 
which  contain  the  actual  surface  geometry,  (2)  bounding  vol¬ 
umes  which  partitions  the  3D  object  space  to  facilitate  ef¬ 
ficient  searching  and  sorting,  and  (3)  surface  relationships 
which  contain  logical  connections  among  surface  and  volume 
sections.  Although  all  three  types  of  nodes  are  applicable 
to  geometrical  modeling  found  in  CAD/CAM  and  computer  graph- 
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ics,  the  last  type  is  also  useful  for  object  matching  and 
recognition  in  computer  vision  applications.  The  system  de¬ 
scribed  here  also  permits  volume  (or  solid)  modeling  as  well 
as  surface  modeling.  Solid  primitives  are  designed  from 
surface  elements  which  enclose  3D  volumes.  Complicated  sol¬ 
id  objects  are  constructed  by  combining  the  solid  primitives 
with  set  operators.  This  allows  the  system  to  be  also  used 
in  constructive  solid  geometry  (CSG)  applications  14,62). 

The  processing  algorithms  which  operate  on  this  surface 
representation  are  based  on  the  method  of  ray-tracing.  This 
method  requires  the  ability  to  traverse  the  data  structure 
and  compute  intersections  between  a  3D  line  (a  ray)  and  3D 
surface  elements.  It  is  an  effective,  versatile,  but  rather 
brute-force  method  which  has  been  used  for  the  generation  of 
shaded  images  (72,64,63)  and  line  drawings  (63),  computa¬ 
tions  of  mass  property  (62,631,  and  conversions  of  represen¬ 
tation  (62,63).  In  the  next  chapter  we  develop  algorithms 
for  (a)  matching  of  3D  surface  descriptions  which  use  ray¬ 
tracing  to  evaluate  a  measure  of  shape  similarity  between 
two  surfaces,  and  (b)  merging  of  overlapping  surface  de¬ 
scriptions  which  use  ray-tracing  to  compute  parameter  trans¬ 
formation,  parameter  projection,  and  conversion  of  represen¬ 
tation. 

In  addition  to  using  the  surface  information  provided 
by  the  photogrammetric  method  described  in  the  previous 
chapter,  surface  data  obtained  from  other  sources  or  syn- 


thetically  designed  may  also  be  entered  into  this  surface 
data  base.  The  section  of  Chapter  3  which  most  directly 
applies  to  the  modeling  process  is  the  conversion  of  the 
surface  information  generated  by  the  reconstruction  method 
into  surface  quadtrees  ”hich  can  then  be  passed  to  the 
various  processing  algorithms.  Chapter  3  also  reviews  other 
forms  of  surface  elements,  in  addition  to  parametric 
patches,  that  are  used  in  the  modeling  system  described 
here,  and  basic  geometric  operations  on  these  elements 
required  by  the  processing  algorithms. 


The  hierarchical  surface  representations  are  stored  in 
directed,  acyclic  graphs.  A  directed  graph  is  generally  de¬ 
fined  to  be  a  finite,  non-empty  set  of  nodes  together  with  a 
set  of  directed  arcs  joining  pairs  of  distinct  initial  and 
final  nodes.  An  acyclic  graph  does  not  contain  any  cycles 
of  arcs.  It  has  one  node  of  in-degree  0  which  is  referred 
to  as  the  root  node.  The  final  nodes  of  arcs  leaving  a  node 
are  called  its  successors :  similarly,  the  initial  nodes  of 
arcs  coming  to  a  node  are  called  its  predecessors.  A  tree 
is  an  acyclic,  directed  graph  in  which  all  nodes  except  the 
root  node  have  in-degree  of  exactly  1.  In  a  hierarchical 
graph  or  tree,  a  node  contains  more  detailed  information 
than  any  of  its  predecessors.  In  an  ordered  graph  or  tree, 
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the  successors  of  a  node  are  listed  in  a  determined  order. 


i 


The  data  structure  of  the  surface  representation,  as 
shown  in  Figures  3.1  and  3.2,  consists  of  the  graph  and  a 
number  of  tables  which  provide  additional  information  about  j 

each  node  or  arc.  The  graph  contains  three  types  of  nodes, 
called  the  R-nodes.  B-nodes.  and  E-nodes,  which  describe 
surface  relations,  bounding  volumes,  and  surface  elements,  j 

respectively.  Although  the  nodes  in  the  graph  are  independ¬ 
ent  of  the  surface  representation,  there  is  a  set  of  sepa¬ 
rate  tables  for  each  type  of  node.  The  actual  information  jj 

about  surface  relations,  bounding  volumes,  and  surface  ele¬ 
ments  is  stored  in  sets  of  R- tables.  B-tables.  and  E-tables, 
respectively.  Detailed  descriptions  of  the  contents  of  I 

these  tables  are  given  in  the  next  three  sections  of  this 
chapter.  The  3D  transformations  represented  by  the  arcs  of 
the  graph  are  stored  in  the  T- table.  A  null  pointer  from  an  | 

arc  implies  identity  transformation.  A  node  entry  in  the 

graph  itself  contains  the  information  shown  in  Figure  3.2, 
and  is  stored  in  a  variable-length  block  in  the  G- table.  I 

This  information  consists  of  (1)  the  type  of  the  node,  (2) 
the  number  of  its  table  and  a  pointer  to  an  entry  in  it,  (3) 
three  flags  which  determine  whether  the  successors  of  the 
node  represent  a  closed  volume  or  an  open  surface,  whether 
they  are  to  be  treated  as  actual  surfaces  and  volumes  or  as 
invisible,  auxiliary  surfaces  and  volumes,  and  whether  they 
contain  a  complete  object  or  a  collection  of  surfaces  and 
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node  type  (B,  E,  or  R) 
table  number 

pointer  to  entry  in  a  table 

"closed  volume  or  open  surface"  flag 

"actual  or  auxiliary  surface"  flag 

"complete  object"  flag 

number  of  successors 

pointer  to  first  successor  (G-table) 

pointer  to  its  transformation  (T-table) 

pointer  to  second  successor  (G-table) 

pointer  to  its  transformation  (T-table) 


pointer  to  last  successor  (G-table) 
pointer  to  its  transformation  (T-table) 


Figure  3.2 


Contents  of  a  graph  node 


volumes,  (4)  the  number  of  successors,  and  (5)  an  ordered 
list  of  pointers  to  the  successors,  each  of  which  has  an 
associated  pointer  to  a  3D  transformation.  Note  that  each 
bounding  volume  or  surface  element  may  potentially  be  de¬ 
fined  in  its  own  local  3D  coordinate  system  while  the  entire 
graph  is  defined  in  a  single  global  object  coordinate  system 
0(x,y,z,w).  Additionally,  the  data  structure  consists  of 
the  A- table  which  contains  surface  attributes  of  individual 
surface  elements,  and  the  I-table  which  contains  identifiers 
to  allow  symbolic  references  to  any  entry  in  any  table.  A 
summary  of  the  graph  representation  is  given  in  Figure  3.3. 

Two  special  cases  of  tree  representation  need  to  be 
mentioned  now:  quadtree  and  octree.  A  quadtree  is  a  hier¬ 
archical,  ordered  tree  in  which  each  non- terminal  node  has 
an  out-degree  of  4.  It  has  been  used  to  encode  efficiently 
information  formatted  on  a  2D  grid,  such  as  digital  images 
[34,62],  and,  therefore,  it  can  also  store  bivariate  surface 
representation  as  will  be  shown  in  Section  3.6.  Similarly, 
an  octree  is  a  hierarchical,  ordered  tree  in  which  each 
non-terminal  node  has  an  out-degree  of  8 .  It  is  being  used 
to  encode  and  process  efficiently  spatial  information  for¬ 
matted  in  a  non-homogeneous  3D  lattice  [41] . 


The  graph 

representation 

employed  here 
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following  main 

advantage  over 

a  tree  representation: 
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allows  a  part 

of  the  surface 
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the  object  space 

to  be  stored  only 

once. 

and 

Node 


Contents 


Typo 


bounding  volume 


surface  element 


surface  relationship 


sphere 
ellipsoid 
parallelepi ped 

plane 

sphere 

quadric  surface 
bicubic  patch 

shape  classification 
alternative  representation 
logical  operation 
object  decomposition 
object  partition 


Tafrls 


transformation 


surface  attributes 
bounding  volumes 
surface  elements 
surface  graph 
symbolic  identifiers 
surface  relationships 
transformations 


Figure  3.3  Summary  of  hierarchical  surface  representation 
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each  of  its  instances  to  be  generated  with  the  appropriate 
3D  transformation  during  a  traversal  of  the  graph.  Simi¬ 
larly,  each  volume  primitive  needs  to  be  defined  only  once 
in  a  standard  orientation  and  then  transformed  into  its  ac¬ 
tual  locations  in  a  solid  object. 

Either  depth-first  or  breadth-first  traversal  of  such  a 
graph  is  possible.  Depth-first  traversal  is  used  in  the 
processing  algorithms  here.  It  requires  a  stack  in  which 
the  pointers  to  all  successors  of  the  current  node,  which 
are  to  be  visited,  are  placed  together  with  their  3D  trans¬ 
formations.  These  transformations  are  computed  by  concate¬ 
nating  the  3D  transformation  of  the  current  node  with  the 
transformation  in  the  arc  pointing  to  the  successor.  The 
next  node  on  the  top  of  the  stack  is  visited  next  until  the 
stack  is  emptied.  However,  if  a  program  frequently  trav¬ 
erses  a  graph,  these  redundant  computations  of  3D  transfor¬ 
mations  can  be  eliminated  by  taking  a  pre-processing  step 
which  (1)  transforms  all  the  coordinate  information  in  the 
graph  to  the  absolute  object  coordinate  system,  and  (2)  con¬ 
verts  the  graph  structure  to  a  tree  structure,  i.e.,  it  gen¬ 
erates  all  instances  of  nodes  with  in-degree  greater  than  1. 

3.2  Relational  Connections 

The  relational  connectivity  of  a  surface  model  gives 
the  logical  and  geometrical  relations  among  the  parts  of  a 
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surface,  among  the  parts  of  an  object/  or  among  a  group  of 
objects.  These  relations  are  mainly  intended  for  any  scene 
analysis  and  object  recognition  work  that  may  be  performed 
with  this  representation.  The  following  five  relational 
categories  are  currently  defined: 

(1)  shape  classification/ 

(2)  alternative  representations/ 

(3)  logical  operations, 

(4)  object  decomposition,  and 

(5)  object  partitions. 

The  relational  connectivity  of  a  surface  is  stored  in  the 
R-nodes  of  the  surface  graph  and  in  the  associated  R-tables. 

The  shape-classification  table  contains  general  infor¬ 
mation  on  the  3D  shape  of  the  first  successor  of  the  current 
R-node.  It  is  classified  into  one  of  these  standard  shapes: 
stick,  blob,  or  box.  The  second  successor  contains  a 
polyhedral  convex-hull  representation. 

The  alter native- representation  table  provides  the  abil¬ 
ity  to  switch  from  one  representation  to  another;  there¬ 
fore,  the  same  surface  may  be  described  by  more  than  one 
representation.  The  first  successor  of  the  current  R-node 
contains  the  original  representation;  the  subsequent  suc¬ 
cessors  contain  the  alternative  representations  which 
usually  are  polyhedral  or  quadric  surface  approximations. 

The  logical- ope rati on  table  defines  logical  relations 
of  the  successors  of  the  current  R-node.  Ac  the  present 
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time,  there  are  defined  three  set  operators:  union.  inter¬ 
section.  and  difference,  which  operate  only  on  closed  vol¬ 
umes.  The  union  operator  produces  the  union  of  all  the  suc¬ 
cessors.  The  intersection  operator  produces  the  inter¬ 
section  of  the  first  successor  with  each  subsequent  succes¬ 
sor.  The  difference  operator  produces  the  difference  be¬ 
tween  the  first  successor  and  each  subsequent  successor. 
Note  also  that  there  is  an  implied  union  of  all  successors 
in  each  B-node.  During  the  ray-tracing  traversal  of  a 
graph,  all  the  intersections  of  the  ray  and  the  volumes  con¬ 
tained  in  the  successor  nodes  are  combined  according  to  the 
specified  set  operator  into  line  segments  of  the  ray  inside 
the  valid  volumes  [63].  These  set  operators  are  used  in  a 
number  of  solid  modeling  systems  [4,62,63]  to  design  compli¬ 
cated  objects  with  CSG  trees  of  volume  primitives.  They  are 
specially  useful  for  interference  checking  and  sectioning. 
In  the  next  chapter,  we  will  propose  a  segmentation  algo¬ 
rithm  that  will  convert  an  arbitrary  3D  solid  object  into  a 
CSG  tree  made  of  these  three  operators  and  a  set  of  solid 
primitives  by  employing  the  surface  matching  algorithm. 

The  obiect-decom do si tion  table  allows  logically  and 
physically  separable  parts  of  an  object  to  be  segmented,  and 
their  spatial  relationships,  such  as  "on  top  of,"  "next  to," 
"supports,"  "attaches,"  "connects,"  "screws  in,"  etc.,  to  be 


identified.  A  list  of  relationships  between  all  pairs  of 
successors  is  provided.  Such  segmentation  is  essential  for 


relational  matching  of  3D  objects  [46,69]. 

The  object- parti tions  table  provides  the  means  to  group 
a  set  of  objects  stored  within  the  same  graph  into  parti- 
tionss  of  similar  objects.  A  successor  of  the  current 
R-node  contains  either  an  object  or  a  group  of  objects  with 
similar  shape.  The  last  successor  of  this  node  then  con¬ 
tains  a  characteristic  description  of  the  objects  in  the 
preceding  successors,  which  typically  is  a  coarse  polygonal 
approximation.  When  a  graph  containing  a  library  of  objects 
is  used  to  identify  an  unknown  object,  the  unknown  object  is 
matched  against  the  objects  in  the  node's  successors  only  if 
the  unknown  object  is  close  to  the  object  description  in  the 
current  node  169] . 

3.3  Hound! gq.  volumes 

A  bounding  volume  is  an  "invisible"  surface  element  or 
a  set  of  surface  elements  completely  enclosing  a  volume 
which  in  turn  may  contain  other  bounding  volumes  or  surface 
elements.  The  purpose  of  bounding  volumes  is  to  partitions 
the  3D  object  space  to  minimize  the  searching  and  sorting  of 
surface  elements  as  is  required  in  a  number  of  algorithms 
for  processing  3D  surfaces.  During  the  ray-tracing  trav¬ 
ersal  of  a  graph,  whenever  a  ray  misses  a  bounding  volume, 
it  also  misses  the  bounding  volumes  and  surface  elements 
within  it  and,  therefore,  the  successors  of  the  bounding 
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volume  are  not  visited.  An  intersection  between  a  ray  and  a 
bounding  volume  is  usually  much  faster  to  compute  than  an 
intersection  between  a  ray  and  a  surface  element.  A  bound¬ 
ing  volume  may  optionally  contain  its  own,  appropriately 
coarse,  description  of  the  surface  elements  within  it, 
usually  computed  by  merging  these  surface  elements.  Should 
a  processing  algorithm  determine,  from,  let  us  say,  the  size 
of  the  bounding  volume  that  a  more  precise  description  of 
the  surface  is  not  required,  it  uses  the  surface  description 
in  the  bounding  volume  and  omits  to  visit  the  bounding  vol¬ 
ume's  successors.  Possible  bounding  volumes  are: 

(1)  sphere, 

(2)  ellipsoid,  and 

(3)  'parallelepiped. 

Ellipsoids  and  parallel pipeds  should  be  oriented  to  minimize 
their  volume.  Only  spherical  bounding  volumes  are  used  in 
the  examples  throughout  this  work. 

A  bounding  volume  is  stored  in  a  B-node  of  the  graph. 
It  contains  the  type  of  the  volume  and  a  pointer  to  an  entry 
in  a  B-table.  A  B- table  contains  the  coefficients  of  the 
volume  and  an  optional  approximation  of  the  surface  elements 
within  the  volume.  There  is  also  a  normal  vector  which  is 
the  average  of  the  normal  vector  field  to  the  surface  ele¬ 
ments  in  the  terminal  nodes  within  the  volume.  The  magni¬ 
tude  of  each  normal  vector  is  proportional  to  the  surface 
area  it  represents.  This,  in  effect,  gives  a  hierarchical 


representation  of  the  surface1 s  Gaussian  Map  [22]. 

A  bounding  sphere  is  computed  by  an  algorithm  for  the 
problem:  "given  a  set  of  N  points  in  3D  space,  find  the 
smallest  sphere  enclosing  them."  The  smallest  enclosing 
sphere  is  defined  either  by  two  points  which  define  its 
diameter  or  by  three  points  of  the  set.  An  algorithm  which 
tests  all  the  spheres  defined  by  two  or  three  points  of  the 
set  and  selects  the  smallest  sphere  which  encloses  all  the 
points  runs  in  0(N2)  time.  An  improved  algorithm  by  Shamos 
[641  locates  the  defining  points  of  the  smallest  sphere  from 
the  extreme  points  of  the  set  by  constructing  a  Voronoi 
diagram  and  runs  in  0 (N  log  N)  time. 

A  bounding  volume  entry  contains  the  following  informa¬ 
tion  stored  in  a  B- table: 

(1)  a  description  of  the  bounding  volume  geometry; 

(2)  an  optional  approximation  to  the  surface  contained 
within  the  volume; 

(3)  error  of  this  approximation,  computed  as  the 
average  distance  between  the  actual  surface  con¬ 
tained  in  the  the  terminal  nodes  and  the  current 
approximation;  and 

(4)  average  surface-normal  vector  integrated  over  the 
actual  surface  within  the  bounding  volume. 

The  error  between  an  actual  surface  and  its  approximation  in 
a  bounding  volume  is  also  computed  by  the  ray-tracing 
method.  Rays  from  the  actual  surface  representation  in  the 
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surface-normal  direction  are  intersected  with  the  approxi¬ 
mated  representation.  The  error  is  computed  as  the  average 
of  the  intersected  distances,  each  weighted  by  the  surface 
area  which  the  ray  represents. 

3.4  Surface  Elements 

A  surface  element  is  a  geometric  representation  of  the 
surface  shape.  Surface  elements  are  contained  in  the 
terminal  E-nodes  of  the  graph  representation.  The  surface 
elements  used  here  are: 

(1)  planar  face, 

(2)  sphere, 

(3)  quadric  surface,  and 

(4)  bicubic  patch. 

The  surface  elements  usually  have  two  different  representa¬ 
tions  -  algebraic  and  parametric.  The  algebraic  representa¬ 
tion  is  an  equation  in  the  form  f(x,y,z,w)  =  0.  This  repre¬ 
sentation,  which  is  used  in  algebraic  geometry,  is  useful  in 
tasks  such  as  "determine  whether  a  given  point  h(x,y,z,w)  is 
on  the  surface,  inside,  or  outside."  The  bounds  of  the  ac¬ 
tual  surface  given  in  this  representation  are  other  algebra¬ 
ic  surfaces  in  structured  logical  trees.  For  a  given  point 
to  be  a  part  of  the  actual  surface,  it  must  satisfy  the 
logical  relations  (inside,  outside)  given  by  the  bounds. 
The  parametric  representation  is  a  vector  g(u,v)  =  [  x(u,v). 
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y(u,v),  z(u,v),  w(u,v)  ]  where  each  component  is  an  equation 
of  the  two  parameters  u  and  v.  This  representation,  which 
is  used  in  differential  geometry  [22],  is  useful  in  tasks 
such  as  "generate  all  points  on  the  given  surface."  The 
bounds  of  a  surface  element  in  this  representation  are 
limits  on  the  parameters  and  y. 

An  E-node  contains  the  type  of  the  surface  element  and 
a  pointer  to  an  entry  in  an  E- table  which  contains  the  coef¬ 
ficients  of  the  surface  representation.  The  two  representa¬ 
tions  are  given  in  detail  for  each  type  of  surface  element 
in  Appendix  C.  There  is  a  separate  E-table  for  each  type  of 
surface  element.  The  surface  bounds  are  defined  in  addi¬ 
tional  auxiliary  tables.  The  planar  elements  have  a  table 
of  edges  and  a  table  of  vertices;  the  quadric  elements  have 
a  table  of  logical  trees  of  other  quadric  surfaces.  Sphere 
representations  do  not  have  bounds,  and  bicubic  patches  are 
defined  over  a  unit  square  of  the  parameter  space.  Each 
surface  element  also  has  a  set  of  surface  attributes  con¬ 
taining  various  properties  of  the  surface  other  than  its 
shape.  These  are  stored  in  the  A- table.  Currently,  only 
reflective,  transmissive,  and  texture  properties  are  mean¬ 
ingful  for  they  can  be  displayed  in  shaded  synthetic  images. 

A  surface  element  in  the  algebraic  representation 
defines  two  half  spaces  -  inside  and  outside.  The  implied 
union  operator  present  in  each  bounding-volume  node  can  be 
used  to  combine  the  inside  half  spaces  into  a  closed  volume 


which  then  can  serve  as  a  volume  primitive.  A  closed  volume 
can  be  defined  with  the  parametric  representation  by  endless 
variations  on  the  isoparametric  brick  theme. 

3.5  Geometric  Operations 

There  are  five  elementary  geometric  operations  that 
need  to  be  performed  on  the  surface  elements  and  the  bound¬ 
ing  volumes: 

(1)  3D  rigid  transformations, 

(2)  evaluation  of  surface-normal  vectors, 

(3)  evaluation  of  surface  curvatures, 

(4)  intersection  with  3D  lines,  and 

(5)  test  for  surface  existence. 

A  3D  rigid  transformation  TCO^-^O^i}  of  a  surface  ele¬ 
ment  or  a  bounding  volume  from  object  space  0^  to  object 
space  0^ i  is  specified  by  a  transformation  matrix  (C.3)  or 
by  a  list  of  x-y-z  displacement,  rotation,  and  scale  parame¬ 
ters  (C.4)  from  which  this  matrix  is  computed  (C.5-8).  An 
entry  in  the  T- table  contains  an  optional  list  of  the  param¬ 
eters  and  the  transformation  matrix  i }  as  well  as 

its  inverse  T{0^ i— >0^} .  The  3D  transformation  is  computed 
during  traversal  of  a  surface  graph  by  concatenating  trans¬ 
formations  stored  in  the  graph's  arcs  until  a  bounding  vol¬ 
ume  or  a  surface  element  is  reached. 

The  normal  vector  to  a  surface  element  or  a  bounding 
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volume  is  computed  from  the  algebraic  representation  by 
equation  (C.10),  or  from  the  parametric  representation  by 
equation  (C.ll).  The  normal-surface  curvature  of  the  param¬ 
etric  representation  is  computed  by  (C.12). 

The  intersection  between  a  surface  element  or  a  bound¬ 
ing  volume  and  a  3D  line,  specified  in  parametric  form  by 
(C.13),  is  computed  by  substituting  the  parametric  line  rep¬ 
resentation  into  the  algebraic  representation  of  the  surface 
element  and  solving  the  general  equation  f(x(t),y(t), 
z(t),w(t))  *  0  for  J;.  Note  that  the  intersection  of  a  line 
and  a  bounding  volume  usually  does  not  have  to  be  solved 
exactly,  it  is  sufficient  to  determine  whether  the  line 
pierces  or  misses  the  bounding  volume.  The  intersections  of 
a  parametic  patch,  which  does  not  have  an  algebraic  repre¬ 
sentation,  and  a  3D  line  is  computed  by  solving  a  system  of 
three  simultaneous  cubic  equations  (C.34)  with  the  parame¬ 
ters  of  the  line  (£) ,  and  the  patch  (u,v)  as  the  unknowns. 
The  algorithm  which  solves  (C.34)  is  given  in  (C.41). 
Alternatively,  this  problem  can  be  presented  as  that  of 
computing  the  intersection  of  a  line  represented  by  the  in¬ 
tersection  of  two  3D  planes  with  the  patch  which  is  a  simi¬ 
larly  complex  problem  of  solving  two  simultaneous  cubic 
equations  with  patch  parameters  (u,v)  as  the  unknowns.  All 
the  results  which  involve  the  computation  of  an  intersection 
between  a  parametric  patch  and  a  line  in  this  report,  either 
in  surface  matching  and  merging  or  in  image  generation,  have 
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been  computed  with  algorithm  (C.41). 

The  intersection  evaluation  is  always  followed  by  a 
check  to  determine  whether  the  intersection  point  is  a  valid 
surface  point,  i.e.,  whether  it  is  within  the  specified  sur¬ 
face  bounds.  For  algebraic  representations  this  involves 
searching  logical  structures  of  bounding  surfaces  also  given 
in  the  algebraic  form.  For  parametric  representations  the 
Ur  £  parameters  of  the  intersection  point  must  be  within  the 
designated  limits  of  the  parameters.  Details  describing  all 
these  geometric  operations  for  each  type  of  surface  element 
are  given  in  Appendix  C. 

3.6  qeneia_tifiiLfl.£. .fi.ucfe.ee.  Quadtrees 


This  section  describes  a  recursive  algorithm  which 
computes  a  quadtree  structure  of  3D  bivariate  surface  data 
formatted  in  a  2D  array.  The  array  is  defined  by  a  2D  pa¬ 
rametric  coordinate  system  P(u,v,w).  We  assume  that  we  are 
given  an  (U+l)  x  (H+l)  array  of  3D  control  points  and, 
optionally,  also  3D  curves  connecting  the  adjacent  points. 
Some  of  the  prints  and  curves  in  the  array  may  be  missing. 
A  surface  patch  can  be  computed  from  a  cycle  of  four  adja¬ 
cent  Points  TiTOf  n(x,y,z,w) ,  hm+lf  n  (x,  y,  z,  w) ,  %+!,  n+1  (x,  y,  z, 
w) ,  and  hjj,,  (x,  y,  z,w) ;  similarly,  it  can  be  computed  from 
a  cycle  of  four  adjacent  boundary  curves.  The  details  of 
these  procedures  are  given  in  Appendix  A.  The  quadtree  al- 


gorithm  uses  the  full  resolution  of  the  array  to  generate 
surface  elements  (pate' es)  and  bounding  volumes  (spheres)  at 
the  bottom  level  of  the  quadtree.  It  then  passes  the  data 
points  at  the  intersections  of  every  other  row  and  column  of 
the  array  to  the  next  level  of  the  quadtree;  that  is,  it 
reduces  the  resolution  of  the  array  by  a  factor  of  two  in 
each  direction.  The  next  level  is  computed  from  this 
reduced  resolution,  linked  to  the  level  below  it,  and  the 
process  is  recursively  repeated  until  the  array  is  reduced 
to  2  x  2  data  points  and  the  root  of  the  tree  is  reached. 

Once  the  coefficients  of  a  patch  are  computed,  a  grid 
of  3D  points  within  the  patch  is  generated.  The  points  are 
used  to  compute  (a)  the  bounding  volume  (a  sphere),  and  (b) 
a  discrete  field  of  normal  vectors  [22],  which  is  averaged 
into  a  single  normal  vector  representing  the  orientation  of 
the  patch.  Each  normal  vector  is  weighted  by  the  surface 
area  it  represents. 

Four  patches  which  share  a  common  control  point  are 
merged  into  a  single  patch  at  the  next  higher  level  of  the 
quadtree,  as  shown  in  Figure  3.4.  The  successors  of  a 
non-terminal  node  in  a  quadtree  are  always  ordered,  relative 
to  the  parametric  coordinate  system  P(u,v,w),  as  is  also 
shown  in  this  figure.  Because  we  need  not  only  the  3D 
positional  data  to  compute  surface  patches,  but  also  the  3D 
slope  data  [Appendix  A],  we  actually  use  the  points  in  a 
reduced  array  to  define  positions,  and  the  deleted  points  to 
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help  in  estimating  slopes.  If  not  all  four  patches  are  de¬ 
fined,  then  only  the  bounding  volume  is  computed  at  the  next 
level.  When  computing  bounding  volumes  at  levels  other  than 
the  bottom  level,  care  must  be  taken  that  each  volume 
encloses  not  only  the  actual  surface  elements  but  also  the 
surface  approximations  in  bounding  volumes  within  the 
current  one. 

An  example  of  the  quadtree  construction  procedure  given 
in  Figure  3.5(a)  shows  a  sheet  of  4  x  4  patches  computed 
from  an  array  of  5  x  5  control  points.  Figure  3.5(b)  shows 
several  bounding  volumes  in  the  quadtree  as  translucent 
spheres.  The  four  smallest  spheres  -  with  green  tint  -  are 
at  the  bottom  level  of  the  quadtree,  each  cast  around  a 
patch.  There  are  12  more  such  spheres,  but  which  are  not 
shown  in  this  image.  The  medium  size  spheres  -  with  blue 
tint  -  are  in  the  middle  level  of  the  quadtree,  each  cast 
around  2x2  patches.  Finally,  the  largest  sphere  -  with 
red  tint  -  is  at  the  top  level  of  the  quadtree,  cast  around 
all  4x4  patches. 

An  outline  of  the  algorithm  which  recursively  builds 
the  quadtree  structure  from  a  sheet  of  surface  data  is  as 
follows : 

(3.1) 

procedure  quadtree (level , M, N) ; 
begin 

if  level  =  0  then 
begin 
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Image  of  a  4  x  4  sheet  of  bicubic  patches  in 
an  object  coordinate  system  0(x,y,z,w):  (a) 
only  surface  elements  shown,  (b)  surface  ele¬ 
ments  and  bounding  volumes  of  a  quadtree  shov*n 


for  m  : =  0  step  1  until  M  do 
begin 

for  n  :=  0  step  1  until  N  do 
begin 

if  surface  data  at  (m,n)  exists  then 
begin 

compute  surface  element  (m,n); 
create  E-node  at  (level, m,n); 
compute  bounding  volume  (m,n); 
create  B-node  at  (level,m,n); 
end? 
end; 
end; 
end; 
else 
begin 

for  m  ;=  0  step  2  mfcil  M  ds 
begin 

for  n  : =  0  step  2  until  N  do 

begin 

if  B-node  defined  at 

(level-1, m, n)  I  (level-1, m+1, n)  I 
(level-1, m,  n+1)  I  (level-1,  m+1,  n+1) 
bben 
begx.n 

merge  bounding  volumes  of  defined  B-nodes 
at  (level-1, m,  n) ,  (level-1, m+1,  n) , 
(level*-l,m,  n+1) ,  and  (level-1, m+1,  n+1) ; 
if  all  4  B-nodes  exist  then 
feegi  n 

compute  surface  approximation; 
compute  error  of  approximation; 
end; 

create  B-node  at  (level , m/2 , n/2 ) ; 
end; 

.and; 

end; 

M  :=  M/2 ; 

N  ;=  N/2; 
end; 

level  ;=  level  +  1; 

if  M>2  I  N>2  then  quadtree (level , M,N) ; 

return; 

end; 


An  alternative  arrangement  to  the  surface  quadtree  rep¬ 
resentation  is  the  binary  tree  representation.  Such  a  tree 
is  created  by  merging  two  adjacent  surface  elements  which 
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share  a  common  boundary  at  any  level  of  the  tree  [Figure 
3.4],  The  direction  of  the  boundary  across  which  the  ele¬ 
ments  are  merged  usually  alternates  from  level  to  level  of 
the  tree.  The  main  goal,  however,  when  selecting  the 
direction  of  a  merger  is  to  minimize  overlap  of  the  bounding 
volumes. 

3.7  Illustrations 

The  example  of  the  reconstructed  test  surface 
'SURFACE. 3'  from  Chapter  1  is  continued  here.  Figure  3.6 
shows  the  quadtree  structure  of  the  surface;  the  branches 
of  the  tree  are  solid  and  the  boundary  curves  of  the  patches 
are  dotted.  The  nodes  of  the  quadtree  are  the  locations  of 


the  centers 

of 

the  bounding 

spheres. 

Figure  3.7  shows 

the 

surface-normal 

vector 

at 

each 

node 

of  the  quadtree; 

the 

branches  of 

the 

tree 

and 

the 

patch 

boundary  curves 

are 

dotted  and  the  normal  vectors  are  solid.  This  model, 
'SURFACE. 3',  consists  of  the  following  information: 

1  sheet 

207  bounding  volumes 
176  control  points 
146  patches  at  level  0 
5  quadtree  levels 

The  surface  shown  in  Figures  3.6  and  3.7  is  contained  in  the 
bottom  level  of  the  quadtree.  Also,  the  four  images  in  Fig¬ 
ure  2.11  were  generated  from  the  bottom  level  of  the  quad¬ 
tree. 
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Quadtree  of  normal  vectors  (normal  vectors 
solid,  tree  branches  and  patches  dotted) 


The  next  three  examples  use  3D  surface  models  which 
were  obtained  from  sources  other  than  our  surface  acquisi¬ 


tion  method.  The  data  for  these  surfaces  had  to  be 
converted  first  into  the  standard  surface  representation. 
The  first  two  examples  demonstrate  objects  modeled  by  large 
numbers  of  sheets  of  bicubic  patches  which  define  closed  3D 
volumes;  the  third  example  represents  open  surface  model 
formed  by  a  single  surface  sheet.  The  first  example  is  a 
surface  description  of  an  F100  engine  blade  designed  by  an 
external  CAD/CAM  system.  The  data  for  the  'F100'  surface 
model  contains: 

47  sheets 

3205  bounding  volumes 
2853  control  points 
2300  patches  at  levels  0 
2  to  7  quadtree  levels  in  a  sheet 

Four  views  of  the  blade  are  shown  in  Figure  3.8,  drawn  as 

line-drawings  of  the  patch  control  points  at  level  1  of  the 

quadtrees.  Four  synthetic  images  of  the  blade,  using  level 

0  of  the  quadtrees,  are  shown  in  Figure  3.9.  In  the  bottom 

two  images  the  blade  is  placed  on  a  planar  surface  and  casts 

a  shadow.  The  second  example  is  a  model  of  a  J79  turbine 

blade.  The  data  for  the  'J79'  model  contains: 

63  sheets 

3819  bounding  volumes 
3551  control  points 
2519  patches  at  level  0 
2-7  quadtree  levels  in  a  sheet 


Four  line  drawings  of  this  blade  are  shown  in  Figure  3.10. 
Four  synthetic  images  of  the  blade,  in  orthographic  pro- 


F100  engine  blade  (3D  network 


Figure  3.9 


Four  images  of  an  F100  engine  blade  (3D  quad¬ 
trees  of  bicubic  patches) 
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jection,  are  shown  in  Figure  3.11.  The  line  drawings  show 
quadtree  data  at  level  1;  the  images  show  quadtree  data  at 
level  0.  The  two  blade  models  shown  in  Figures  3.9  and  3.11 
were  assigned  highly  specular  reflections  to  simulate  a 
metalic  surface. 

The  third  example  represents  a  3D  terrain  model  ob¬ 
tained  from  a  digital  contour  map  of  the  Fall  River  Pass 
Quadrangle  in  Colorado.  The  actual  area  of  the  modeled 
terrain  is  6.6  x  8.5  miles.  The  ratio  of  the  vertical  scale 
to  the  horizontal  scales  in  this  model  is  7:1:1  which 
significantly  exaggerates  the  terrain's  elevation.  The  sur¬ 
face  model  'TERRAIN'  contains  the  following  data: 

1  sheet 

5489  bounding  volumes 
4218  control  points 
4088  patches  at  level  0 
8  quadtree  levels 

There  are  four  views  of  the  terrain,  drawn  at  level  2  of  the 
quadtree,  shown  in  Figure  3.12.  There  are  four  similar 
views  of  the  terrain  shown  in  the  images  of  Figure  3.13. 
The  first  image  (upper  left)  was  generated  from  data  at 
level  0  of  the  quadtree  which  contains  56  x  76  patches.  An 
image  of  the  original  contour  map  is  shown  under  the 
terrain,  mapped  on  a  rectangular  plane.  The  other  three  im¬ 
ages  were  generated  from  data  at  level  2  of  the  quadtree 
which  contains  14  x  19  patches.  Notice  the  difference  in 
detail  between  the  two  representations.  The  terrain  is 
shown  from  SSW  (south-southwest),  SSE,  NNE,  and  NNW,  respec- 
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Figure  3.13  Four  images  of  a  terrain  model  (3D  quadtree  of 
bicubic  patches) 
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tively,  with  the  sun  at  SE,  30  degrees  above  the  horizon. 
The  colors  of  the  surface  are  computed  as  a  function  of 
elevation  from  a  look-up  table;  they  range  from  dark  green 
in  the  valleys  through  light  green,  yellow,  orange,  and 
light  brown  to  dark  brown  at  the  mountain  peaks. 

The  solid  modeling  capabilities  of  this  system  are 
demonstrated  in  the  last  two  examples  given  in  this  chapter. 
First,  the  logical  set  operations  are  shown  on  two  volume 
primitives  [Figure  3.14] .  Volume  primitive  'A'  is  formed  by 
union  of  three  surface  half-spaces:  two  spheres  and  one 
cylinder;  similarly,  primitive  1 B'  is  formed  by  union  of 
two  planes  and  an  ellipsoid  [Figure  3.14(a)].  The  union, 
difference,  and  intersection  operations,  applied  to  the  two 
volume  primitives,  are  shown  in  Figures  3.14(b),  (c),  and 
(d),  respectively. 

Finally,  two  section  views  of  the  J79  blade  are  shown 
in  Figure  3.15.  They  are  made  by  intersection  and  differ¬ 
ence,  respectively,  of  the  blade  model  and  an  invisible  par¬ 
allelepiped  which  encloses  the  left  half  of  the  blade.  One 
face  of  the  parallelepiped  intersects  the  blade  in  the 
vertical  and  front-back  directions.  Since  the  parallelepi¬ 
ped  is  invisible  the  sectioned  blade  appears  hollow;  the 
very  dark  surface  is  the  back  surface  of  the  blade  in  shadow 
of  the  front  surface. 
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3.8  Summary 


In  this  chapter  we  presented  an  experimental  modeling 
system  for  hierarchical  representation  of  3D  surfaces,  vol¬ 
umes,  and  objects.  The  system  is  intended  for  a  wide  range 
of  applications,  from  CAD/CAM  to  computer  vision.  Desirable 
features  such  as  good  user  interface  for  designing  new  sur¬ 
faces  and  editing  existing  surfaces,  or  automatic  testing 
for  validity  and  consistency  of  representation  (detection  of 
nonsense  objects)  that  are  mandatory  in  an  actual  modeling 
system  are  beyond  the  scope  of  this  work.  All  application 
algorithms  which  process  this  surface  representation  use  a 
common  ray-tracing  traversal  procedure  to  obtain  a  list  of 
intersections  of  a  ray  with  the  surfaces  in  a  directed, 
acyclic  graph.  Surface  data  is  usually  provided  to  the  sys¬ 
tem  as  bivariate  surface  elements  which  are  then  structured 
into  hierarchical  quadtrees.  However,  representation  by 
polyhedral  and  quadric  surfaces  is  also  available.  New  sur¬ 
face  elements  can  be  simply  added  to  the  system,  provided 
that  the  five  geometric  operations  (Section  3.5)  are  de¬ 
fined.  These  are  the  only  operations  in  the  entire  system 
which  depend  on  the  type  of  surface  representation. 
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CHAPTER  4 


SURFACE  MATCHING 

Following  the  prerequisite  processing  described  in  the 
last  two  chapters,  we  finally  reach  the  heart  of  this  work  - 
matching  of  3D  surfaces.  In  the  context  of  this  work, 
matching  of  3D  surfaces  refers  exclusively  to  finding  a 
spatial  registration  of  two  3D  surface  descriptions  that 
maximizes  their  shape  similarities.  If  a  measurement  of  the 
shape  similarities  is  acceptable  then  the  surfaces  are  said 
to  match,  and  the  spatial  registration  aligns  them.  In  or¬ 
der  to  find  such  a  match,  the  two  surface  descriptions, must, 
of  course,  completely  or  partially  overlap,  i.e.,  there  must 
be  a  set  of  surface  points  common  to  both  surface  descrip¬ 
tions.  Note  that  we  do  not  refer  here  to  other,  possibly 
similar,  types  of  matching  such  as  (a)  matching  2D  projec¬ 
tions  with  2D  models  of  3D  objects  [17,29,14],  (b)  matching 
2D  projections  with  3D  models  [9,11],  or  (c)  relational 
matching  of  features  of  3D  objects  and  3D  models  [69]. 

There  are  a  number  of  problems  where  two  or  more  3D 
surfaces  need  to  be  matched.  For  example,  a  partial  3D  sur¬ 
face  description  of  an  unknown  object  is  obtained  from  one 
(stereo)  vantage  point.  The  object  is  to  be  recognized  by 
matching  this  partial  3D  surface  description  with  a  set  of 
3D  models  of  various  objects.  Or  pieces  of  broken  pottery 
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or  fragments  of  objects  need  to  be  matched  to  allow  recon¬ 
struction  of  the  original  objects.  Or  a  3D  description  of 
surface  below  an  aircraft  needs  to  be  matched  with  a  3D 
terrain  model  to  determine  accurately  the  position  of  the 
aircraft.  In  the  scope  of  this  work,  we  need  to  match  3D 
surface  segments  of  an  object,  which  were  generated  while 
the  object  was  in  several  stable  positions  and  parametrized 
by  several  parameter  systems,  in  order  to  build  a  complete 
3D  model  of  the  object. 

The  matching  algorithm  of  two  3D  surfaces,  which  may 
completely  or  partially  overlap,  consists  of  two  major 
steps:  (1)  initial  estimates  of  the  surface  registration 
are  computed  by  alignment  of  known  points  on  both  surfaces 
or  alignment  of  surface-normal  vectors  representing  surface 
orientations,  and  (2)  a  heuristic  search  improves  these  es¬ 
timates  by  varying  transformation  parameters  to  find  an 
acceptable  solution.  The  measurement  of  shape  similarity 
between  the  two  surfaces  is  computed  by  an  evaluation  func¬ 
tion  from  the  following  information,  obtained  for  a  number 
of  points  distributed  on  both  surfaces:  (1)  Euclidean 
distance  to  the  other  surface,  (2)  angular  difference  be¬ 
tween  normal  vectors,  and  (3)  difference  in  local  surface 
curvatures.  The  ray-tracing  traversal  of  a  surface  descrip¬ 
tion  computes  this  information.  Given  a  point  on  one  sur¬ 
face  and  the  surface-normal  vector  at  this  point  it  finds 
(1)  the  nearest  intersection  with  the  other  surface,  (2)  the 
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surface-normal  vector  at  the  point  of  intersection,  and  (3) 
the  surface  curvature  at  the  point  of  intersection.  The 
matching  algorithm  is  independent  of  the  surface  representa¬ 
tion  which  is  confined  only  to  the  ray-tracing  process. 

Following  the  description  of  the  matching  problem,  we 
present  a  number  of  algorithms  which  merge,  into  a  single 
concise  description,  those  overlapping  surfaces  that  were 
either  obtained  in  the  same  object  coordinate  system,  or 
transformed  into  the  same  object  coordinate  system  by  the 
matching  algorithm.  The  first  algorithm  merges  two  param¬ 
etrized  surfaces  into  a  common  2D  parameter  coordinate  sys¬ 
tem  by  a  transformation  of  one  parameter  coordinate  system. 
The  second  algorithm  and  its  variations  reparametrize  the 
surface  descriptions  into  a  new  parameter  coordinate  system 
by  projecting  this  system  on  the  surfaces.  They  can  also  be 
used  to  parametrize  surfaces  given  only  in  algebraic  repre¬ 
sentations.  The  third  and  last  merging  algorithm,  which  is 
applicable  only  to  surface  descriptions  of  a  solid  object, 
converts  the  surface  representations  into  a  polyhedral  vol¬ 
ume  representation  made  of  rectangular  parallelepipeds. 
This  solid  representation  can  be  further  converted  into  the 
octree  representation  144]. 

Finally,  this  chapter  closes  with  descriptions  of  three 
applications  of  the  matching  algorithm  to  (1)  generation  of 
complete  3D  models,  (2)  surface  and  object  recognition,  and 
(3)  surface  and  volume  segmentation  with  surface  and  volume 


primitives.  The  first  application  also  uses  the  merging  al¬ 
gorithms.  It  has  been  implemented  and  tested  with  several 
objects.  The  last  two  applications  are  described  as 
suggested  approaches  to  these  problems. 


There  are  two  basic  types  of  surface  matching  [Figure 
4.1]  that  are  of  interest  in  this  work.  Given  two  3D  sur¬ 
face  descriptions,  represented  by  graphs  £  and  £' : 

(a)  surface  description  £  is  completely  contained  in 

description  £'  [Figure  4.1(a)];  that  is,  the  in¬ 
tersection  set  of  £  and  £'  is  £  (e.g.,  a  surface 

segment  is  being  matched  with  a  complete  object 
model ) ;  and 

(b)  surface  description  £  is  only  partially  contained 
in  description  £'  [Figure  4.1(b)];  that  is,  the 
intersection  set  of  £  and  £'  is  a  subset  of  £  and  a 
subset  of  £'  (e.g.,  two  overlapping  surface  seg¬ 
ments  of  an  object  are  being  matched). 

These  two  types  of  matching  differ  in  the  approach  that  the 
matching  algorithm  uses  in  generating  the  initial  estimates 
of  the  surface  registrations.  In  case  (a),  it  emphasizes 
the  spatial  registrations  that  completely  match  £  to  £'.  In 
case  (b),  it  emphasizes  the  spatial  registrations  that  match 
only  a  portion  of  £  to  a  portion  of  £'.  The  first  type  of 
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■B^!l 


Figure  4.1 


Types  of  surface  matching:  (a)  5  completely 

overlaps  (b)  £  partially  overlaps  £' 


surface  matching  is  used  for  object  recognition.  The  second 
type  is  used  for  generation  of  complete  object  models,  and 
for  their  segmentation  with  surface  and  volume  primitives. 

A  similar  type  of  surface  matching  could  also  be  used 
for  solving  the  3D  "jig-saw  puzzle"  problem  where  surface 
segments  with  only  common  boundary  curves  (e.g.r  broken 
pottery  pieces)  need  to  be  matched  and  spatially  registered. 
In  this  problem  each  surface  segment  would  be  extrapolated 
and  labeled  as  either  "actual"  or  "extrapolated"  sections. 
The  matching  algorithm  would  then  compute  spatial  regis¬ 
trations  by  matching  the  "actual"  sections  of  G  with  the 
"extrapolated"  sections  of  £',  and  the  "extrapolated"  sec¬ 
tions  of  Q  with  the  "actual"  sections  of  £'. 

4.2  Surface  Transformations 

The  goal  of  the  matching  algorithm  is  to  compute  a  ri¬ 
gid  3D  transformation  that  matches  two  surface  descriptions 
and  aligns  them  into  the  same  object  coordinate  system.  In 
this  section  we  describe  two  types  of  rigid  3D  transforma¬ 
tions  [Figure  4.2J,  used  by  the  matching  algorithm,  and  the 
computational  methods  that  generate  them.  The  two  transfor¬ 
mations  consist  of  different  numbers  of  transformation  pa¬ 
rameters.  The  first  transformation  is  a  general  transforma¬ 
tion  between  two  arbitrarily-oriented  object  coordinate  sys¬ 
tems.  It  consists  of  three  translation,  three  rotation, 


88 


and,  optionally,  three  scale  parameters.  The  second  trans¬ 
formation  is  a  limited  transformation  between  two  object  co¬ 
ordinate  systems  which  have  one  of  the  three  axes  parallel. 
It  consists  of  two  translation,  one  rotation,  and, 
optionally,  two  scale  parameters.  It  is  used  only  when 
there  is  k  priori  knowledge  that  the  object,  whose  surface 
segments  are  being  matched,  remained  in  the  same  stable 
position.  The  scale  parameters  are  normally  not  used  since 
it  is  assumed  that  all  the  surface  descriptions  have  been 
obtained  to  the  same  scale  -  the  actual  size  of  the  sur¬ 
faces.  They  would  only  be  used  for  object  segmentation  with 
surface  and  volume  primitives  defined  at  some  standard  size. 

A  general  transformation  [Figure  4.2(a)] 
from  homogeneous  object  coordinate  system  0(x,y,z,w)k  to  ho¬ 
mogeneous  object  coordinate  system  0(x,y,z,w)ki  is  expressed 
as: 

(4.1) 


This  transformation  consists  of  the  following  transformation 
parameters:  three  translations  (tx,  ty,  tz),  three  rota¬ 
tions  (rx,  ry,  rz),  and,  optionally,  three  scale  factors 
(sx,  sy,  sz) .  The  twelve  unknown  coefficients  of  T(Oi.— « ) 


can  be  determined  from  the  coordinates  of  four  non-coplanar 
points  in  0(x,y,z,w)k  and  their  corresponding  images  in 
0(x,y,z,w)ki . 

A  limited  3D  transformation  T(0k— »0kt}  [Figure  4.2(b)] 
from  homogeneous  object  coordinate  system  0(x,y,z,w)k  to  ho¬ 
mogeneous  object  coordinate  system  0 (x, y, z, w) k ,  with  zk  = 
zKi  is  expressed  as: 

(4. 
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This  transformation  consists  of  these  transformation  parame¬ 
ters:  two  translations  (tx,  ty) ,  one  rotation  (rz)r  and, 
optionally,  two  scale  factors  (sx,  sy) .  The  six  unknown  co¬ 
efficients  of  this  T{0k— »0ki)  can  be  determined  from  the  js 
and  y  coordinates  of  three  non-colinear  points  in 
0(x,y,z,w)k  and  their  corresponding  images  in  0  (x,y, z,w) k 1 . 
This  transformation  is  available  to  simplify  the  matching 
process  when  the  two  surface  segments  being  matched  have 
been  obtained  from  an  object  in  the  same  stable  position 
with  respect  to  the  x-y  planes  in  the  0(x,y,z,w)k  and 
0(x,y,z,w)k«  coordinate  systems. 

In  general,  a  rigid  3D  transformation,  composed  of 
translation  and  rotation,  is  defined  as: 
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or,  for  short, 

Ti  »  R  "H  +  t  (4.4) 
where  t  is  a  translation  (displacement)  transformation  and  R 
is  an  orthogonal  rotation  transformation  which  imposes  the 
condition: 

R  e  •  r!  =  e  •  f  (4.5) 
for  all  vectors  e  and  1  in  0(x,y,z,w).  This  condition 
preserves  the  distance  between  any  two  points  and  the  angle 
between  any  two  vectors.  In  addition,  the  determinant  of  R 
must  be  positive  to  preserve  orientation  (i.e.  to  avoid  re¬ 
flections)  in  the  right-handed  object  coordinate  system 
0(x,y,z,w).  The  inverse,  R~*,  of  an  orthogonal  matrix  R  is 
equal  to  its  transpose  Rt.  Note  that  R  is  the  upper-left  3 
x  3  principal  submatrix  of  T  in  equation  (4.1),  and  t  is  the 
upper-right  3x1  submatrix  of  T.  A  transformation  which 
contains  scale  factors  other  than  unity  is  not  a  true  rigid 
transformation  because  it  does  not  preserve  distances  but 
only  angles. 

The  matching  algorithm  uses  two  methods  to  compute  the 
object  transformation  matrix  (4.1)  (the  limited  transforma¬ 
tion  (4.2)  is  only  a  special  case  of  (4.1)).  The  first 
method  composes  a  transformation  matrix  for  each  transforma- 
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tion  parameter  and  multiplies  these  matrices  into  the  total 
transformation  [Appendix  C.1.21.  A  rotation  matrix  composed 
by  this  method  is  guaranteed  to  be  always  orthogonal. 


The  second  method  determines  the  transformation  from 
matched  points  in  the  two  object  coordinate  systems.  Given 
a  point  Ti(x,y,  z,w)  in  0(x,y,z,w)k  and  its  image  h'(x,y,z,w) 
in  0(x,y,z,w)|ci  we  rewrite  (4.1)  as  a  system  of  three  linear 
equations  with  twelve  unknowns: 

(4.6) 
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or,  for  short, 

H  (4.7) 
With  three  additional  points  we  can  express  system  (4.6)  as 
twelve  equations  with  twelve  unknown  coefficients  of  T.  In 
general,  for  N  2  4  points  the  least-squares  method  is  used 
to  minimize  the  sum  of  squares  of  distances  between  the 
matched  pairs  of  points  'Hi  and  h'^  in  equation  (4.4): 


N 

minC£  |  H'i  -  R  T»i  +  t  |2]  (4.8) 

i=l 


which  yields  the  standard  least-squares  solution: 


This  solution,  however,  does  not  guarantee  R  to  be  orthogo¬ 
nal.  Once  we  have  obtained  the  translation  vector  7  from 
(4.9),  we  can  improve  the  othogonality  of  R  by  also  using 
the  inverse  transformation  of  (4.4),  namely: 

H  ■  R"1  (K'-t)  =  RtrtJ'-I)  (4.10) 

to  minimize  the  squares  of  distances,  as  in  (4.8),  with  the 

transformation  R  and  its  inverse  Rfc: 

N 

min£^  |  h'i-R  Tij  +  t  |2  +  \  Rfc  (h1  j-t) -T^  |2^]  (4.11) 

i*l 

The  solution  to  (4.11)  is  similar  to  (4.9),  and  still  not 
necessarily  orthogonal.  Note  that  to  obtain  an  exactly  or¬ 
thogonal  R  we  would  have  to  solve  a  system  of  nine  simulta¬ 
neous  quadratic  equations.  To  avoid  this,  the  matrix  R  ob¬ 
tained  from  (4.11)  is  converted  to  an  orthogonal  matrix  with 
the  following  procedure: 

(4.12) 

(1)  convert  each  row  (column)  vector,  7^,  of  R  to  a 
unit  vector; 

(2)  find  two  most  orthogonal  row  (column)  vectors,  7 i 


(3)  make  r^  and  rj  orthogonal: 

■  *i  -  (Ti-tj)  It  j ; 

(4)  convert  new  7^  to  a  unit  vector; 

(5)  make  7^  othogonal  to  7A  and  7j: 

7k  -  ?i  X  tjj 

(6)  preserve  orientation: 

if  | R |  <  0  then  7^  =  -7^. 

In  summaryf  a  rigid  3D  transformation  is  computed  from 
the  coordinates  of  points  in  0(x,yrz,w)  and  their  images 
in  O' (x,y,z,w)  by  equation  (4.8)  which  obtains  translation 
t,  equation  (4.11)  which  obtains  rotation  R,  and  procedure 
(4.12)  which  makes  R  orthogonal.  The  last  step  prevents 
skewed  and  sheared  transformations. 

4.3  Matching  Algorithm 

The  matching  algorithm  computes  a  3D  rigid  transforma¬ 
tion  T,  composed  of  a  rotation  matrix  R  and  a  translation 
vector  t,  which  spatially  registers  (aligns)  surface  de¬ 
scription  £  with  surface  description  £•.  This  transforma¬ 
tion  occurs  when  the  orientation  and  shape  differences  be¬ 
tween  the  two  surfaces  are  minimized  (i.e.,  the  shape  simi¬ 
larity  is  maximized).  The  orientation  and  shape  differences 
are  evaluated  at  a  number  of  evaluation  points  on  surface  £ 
from  information  provided  by  the  ray-tracing  procedure 
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[Chapter  3].  For  a  ray  in  the  direction  of  the  sur¬ 
face-normal  vector  at  an  evaluation  point,  Ti(x,y,z,w)^r  on 
surface  £  this  procedure  finds  the  nearest  intersection  with 
fi',  the  surface-normal  vector,  and  the  normal  surface  curva¬ 
ture,  both  evaluated  at  the  point  of  intersection.  The  dif¬ 
ferences  in  surface  orientation  and  shape  are  computed  from 
these  components: 

(1)  position  difference.  (T) ,  is  the  3D  Euclidean 
distance  in  the  object  coordinate  system  O' (x,y, z,w) : 

p^  (T)  =  |  T  'E(x,y,z,w)i  -  h'  (x,y,z,w>i  |  (4.13) 

where 

T  -  current  transformation  from  0(x,y,z,w) 

to  O' (x,y,z,w) 

h(x,y,z,w)i  *  an  evaluation  point  on  surface  fi 

"H •  (x,y,z,w>i  =  the  point  of  intersection  on  surface  £' 

(2)  orientation  difference,  a^  (T) ,  is  the  angular  dif¬ 
ference  of  the  surface  unit-normal  vectors  oriented  to  the 
"outside"  of  the  surfaces: 

at  (T)  «  |  (R  Nj)  •  N'i  -  1.0  1  (4.14) 

where 

R  *  current  rotation  from  0(x,y,z,w)  to  O'  (x,y,z,w) 

Ni  «  unit  normal  vector  at  h(x,y,z,w)i 

N'i  =  unit  normal  vector  at  "H '  (x, y, z,w) ^ 
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(3)  curvature  difference*  c^ (T) ,  is  the  magnitude  of 
difference  of  the  surface  normal  curvatures: 

C£  (T)  «  I  <3i  ~  q'i  |  (4.15) 

where 

■  normal  surface  curvature  at  h(x,y,  z,w)^ 
q'i  «  normal  surface  curvature  at  "H*  (x, y,z,w)  ^ 

The  total  surface  difference,  d^ (T) ,  at  this  point  is 
then  defined  as  a  linear  combination  of  these  three 
components : 

d^ (T)  *  WpPi (T)  +  waai (T)  +  wcCi  (T)  (4.16) 

where 

wp  ■  a  weight  of  position  distance  p^ (T) 
wa  *  a  weight  of  angular  difference  a^ (T) 
wc  =  a  weight  of  curvature  difference  Ci (T) 

Given  N  evaluation  points  on  surface  G,  which  are  used 
in  computing  the  surface  registration,  we  seek  a  3D  trans¬ 
formation  7  such  that: 

N 

D  (f )  w^di  (T)  <  epsilon  (4.17) 

i-1 

where 

d|  (T)  *  surface  difference  at  point  'E(x,y,z,w)i 

w^  -a  weight  of  point  h(x,y,z,w)^ 

epsilon  ■  an  acceptable  difference  between  the  two  sur¬ 
face  descriptions  G  and  £'  for  a  match 
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An  additional  condition  requires  that  a  correct  surface 
match  be  found  for  at  least  U  of  the  H  evaluation  points 
where  4  «  JJ  i  Notice  that  a  ray  from  Ti(x,y,z,w)  i  may 
(a)  miss  in  which  case  d^  =  0,  or  (b)  find  an  incorrect 

match  with  £'  (Figure  4.3],  in  which  case: 

(4.18) 

|h' (x,y, z,w) j-h' (x,y, z,w) j |  >>  | h (x, y, z, w) ^-h (x, y, z, w) j | 
rather  than 

(4.19) 

|H' (x,yf ZrwJi-h' (x,y,z,w) j |  •  |h(x,y, z,w)i-h(x,y,z,w) j | 

for  any  point  Ti (x, y,z,w)  j,  near  li(x,y,z,w)  and  its  inter¬ 
section  "H '  (x,y,z,w)  j  on  £'  (i.e.,  as  stated  in  Section  4.2: 
the  distance  of  any  two  points  must  be  preserved  by  T) .  If 
the  relationship  (4.18)  is  valid  for  the  point  H(x,y,z,w)^ 
then  d^  «  0. 

Having  defined  the  method  which  evaluates  the  surface 
difference,  we  shall  now  describe  the  actual  search  algo¬ 
rithm  for  a  spatial  registration  which  minimizes  this  dif¬ 
ference  between  the  two  surfaces.  A  basic  matching  algo¬ 
rithm  which  would  blindly  wander  in  its  task  without  any 
guidance  can  be  outlined  as: 

(4.20) 

A  naive  algorithm 


(1)  Select  a  set  of  H  evaluation  points  on  £. 


(2)  Generate  a  new  transformation  T  of  £. 

(3)  Compute  surface  difference,  d^,  at  each  evaluation 
point,  h(x,y,z,w)^,  transformed  by  T,  with  £'. 

(4)  If  the  total  difference  D  (T)  is  less  than  epsilon 
and  g  is  much  greater  than  4  then  the  surfaces 
match  and  T  aligns  them,  otherwise  go  to  (2). 

Clearly,  an  exhaustive  search  for  the  transformation 
parameters  of  T  would  be  computationally  prohibitive.  A 
hill-climbing  method,  that  would  adjust  each  parameter  in 
turn,  would  also  be  computationally  excessive  and  probably 
would  fail  in  a  number  of  cases.  A  viable  method  to  reduce 
the  search  task  'is  to  use  a  state-space  search  method  [48, 
49]  with  an  evaluation  function  to  guide  the  search  for  the 
transformation  parameters  of  T  that  satisfy  (4.17).  This 
method  represents  the  search  process  by  a  graph  whose  nodes 
are  generated  by  successor  operators  which  attempt  to 
improve  the  node's  surface  registration.  An  evaluation 
function,  £,  provides  a  ranking  of  the  graph  nodes  to 
determine  which  nodes  are  most  likely  to  be  on  a  path  to  the 
solution  and  should,  therefore,  be  expanded.  An  outline  of 
this  algorithm,  adapted  to  the  domain  of  our  problem,  can  be 
specified  as  follows: 

(4.21) 

An  ordered-search  algorithm 

(1)  Generate  start  nodes,  and  their  T  transforma^- 


tions;  put  them  on  list  OPEN ,  compute  e(n^)  for 
each  start  node  ji^. 

(2)  If  OPEN  is  empty  exit  with  failure. 

(3)  Remove  from  OPEN  a  node  n  with  the  smallest  e(n) 
and  put  it  on  list  CLOSED.  If  CLOSED  becomes  full 
exit  with  failure. 

(4)  If  D  (T)  of  n  is  less  than  epsilon  and  M  of  jj  is 
much  greater  than  4,  exit  with  solution  T  of  jj. 

(5)  Expand  node  n,  generate  all  of  its  successors  by 
computing  their  T  transformations;  for  each  suc¬ 
cessor,  nif  compute  eto^. 

(6)  Put  all  successors,  which  are  not  already  on  OPEN 
or  CLOSED,  on  OPEN  and  link  them  to  n.  If  OPEN 
becomes  full  exit  with  failure. 

(7)  For  all  successors  which  are  already  on  OPEN  or 
CLOSED,  if  new  e(n^)  is  lower  than  the  old  etn^ 
then  replace  it;  move  from  CLOSED  to  OPEN  all 
nodes  whose  value  eCn^  was  lowered.  If  OPEN 
becomes  full  exit  with  failure. 

(8)  Go  to  (2) . 

Each  start  node  is  the  root  node  of  a  search  tree;  all  the 
terminal  nodes  of  a  search  tree  are  on  list  OPEN.  The 
critical  parts  of  this  algorithm  are  the  choice  of  the  eval¬ 
uation  function  e(n),  the  test  for  the  presence  of  a  succes¬ 
sor  node  on  OPEN  or  CLOSED  in  steps  (6)  and  (7),  the  selec- 


tion  of  the  evaluation  points  on  surface  G  to  be  matched 
with  £',  and  the  strategies  of  the  successor  operators  which 
(a)  generate  start  (root)  nodes  of  the  search  trees  in  step 
(1),  and  (b)  expand  (terminal)  nodes  on  OPEN  in  step  (5). 
These  parts  of  the  matching  algorithm  shall  now  be  described 
in  detail. 

The  evaluation  function  critically  affects  the  search 
process.  A  function  which  is  too  generous  will  cause  the 
expansion  of  too  many  nodes.  On  the  other  hand,  a  function 
which  ignores  the  potential  of  some  nodes  can  lead  to  a  fu¬ 
tile  search.  The  currently  used  evaluation  function,  e(n), 
is  a  weighted  sum  of  the  surface  difference  D (T) ,  the 
inverse  of  the  number  of  matched  evaluation  points,  and  the 
length  of  the  tree  path  from  the  current  node  to  its  start 
node : 

e  (n)  =  wdD(T)  +  wm/M  +  wLL(n)  (4.22) 

where 

D  (T)  *  difference  of  £,  transformed  by  T,  and  £' 
wD  =*  a  weight  of  surface  difference  D  (T) 

M  ■  number  of  matched  evaluation  points 

wM  *  a  weight  of  inverse  of  number  of  matched  points 

L(n)  *  the  path  lenght  from  jj  to  its  start  node 
wL  «  a  weight  of  path  lenght  L(n) 

The  first  two  terms  of  the  evaluation  function  (4.22) 
provide  heuristic  information  about  the  quality  of  the 


spatial  registration 


The  ordered-search  algorithm  has  to  check,  in  steps  (6) 
and  (7),  whether  a  successor  node  is  already  present  on  OPEN 
or  CLOSED.  Each  node  contains  the  coordinates  of  a  standard 
3D  vector,  H,  transformed  by  the  node's  T  transformation. 
If  the  distance  of  the  transformed  vector,  T  I,  of  a  succes¬ 
sor  node  is  less  than  a  selected  tolerance  from  T  H  of  a 
node  already  on  OPEN  or  CLOSED,  then  the  transformation  of 
the  successor  node  is  assumed  to  be  the  same  as  that  of  the 
node  on  OPEN  or  CLOSED  and  the  successor  node  is  already 
present  there. 

The  set  of  evaluation  points  on  surface  G,  where  the 
surface  difference  (4.17)  is  computed,  should  fairly  repre¬ 
sent  the  shape  of  the  surface.  There  are  two  possible 
approaches,  considered  here,  to  the  selection  of  these  eval¬ 
uation  points:  (a)  find  "critical"  points  at  high  surface 
curvature  (e.g.,  surface  vertices  and  edges),  or  (b)  gener¬ 
ate  regularly  spaced  points  on  a  parametric  grid.  The  first 
approach,  although  potentially  more  powerful,  suffers  from  a 
number  of  problems.  It  is  difficult  to  find  individual  sur¬ 
face  points  with  significantly  higher  curvature  on 
relatively  smooth,  multiply-curved  surface  segments  (see 
'SURFACE. 3'  in  Figures  2.9-11).  On  surfaces  with  constant 
curvature  (e.g.,  spheres  and  cylinders)  this  method  fails 
completely.  On  highly-curved  surfaces  these  points  tend  to 
cluster  along  surface  edges  and  be  colinear.  A  search  to 


[•] 


3 


locate  points  with  maximum  surface  curvature  is  also 
computationally  demanding. 

In  the  second  approach,  adapted  by  the  matching  algo¬ 
rithm,  the  evaluation  points  are  regularly  spaced  on  a  sur¬ 
face  and  hierarchically  structured.  A  natural  representa¬ 
tion  for  this  approach  is  the  surface  quadtree  developed  in 
Section  3.6.  The  control  (corner)  points  of  parametric 
patches  are  used  as  the  evaluation  points.  In  the  initial 
stages  of  a  search,  the  algorithm  uses  the  control  points 
from  the  high  levels  of  a  quadtree  £;  as  the  value  of  D (T) 
decreases,  the  control  points  are  replaced  by  points  from 
lower  levels  of  the  quadtree  and  their  number,  therefore, 
increases.  Similarly,  initially  in  a  search  the  algorithm 

t 

intersects  rays  from  these  evaluation  points  with  patches 
stored  in  the  high  levels  of  quadtree  £'.  As  the  value  of 
D (T)  descreases,  patches  from  lower  levels  are  intersected. 
In  summary,  the  matching  algorithm,  therefore,  uses  the  hi¬ 
erarchical  surface  description  as  follows:  while  the  esti¬ 
mate  of  the  spatial  registration  is  coarse,  only  a  coarse 
surface  description  is  used  to  evaluate  D  (T) ;  as  the 
spatial  registration  improves,  D  (T)  is  evaluated  at  more 
points  with  more  accurate  surface  description.  Selection  of 
the  quadtree  level  used  in  the  evaluation  of  the  surface 
difference  is  a  function  of  D(T): 


0 


I  if  DC!)  i  ‘W1* 

level  i  if  Dmin(i)  <  D  (T )  <  Dmax(i>  (4.23) 

levelmax  if  DmaK!1*«lmax-1>  <  D  (T) 

where  0  <  i  <  level^^,  and  levelmax  is  the  highest  level  of 
a  quadtree  containing  a  surface  approximation.  An  evalu¬ 
ation  point  may  be  assigned  a  weight  proportional  to  its 
surface  curvature.  This  weight  is  used  in  (a)  evaluation  of 
D(T),  and  (b)  generation  of  new  T  when  a  graph  node  is  being 
expanded.  For  surface  and  volume  primitives  it  may  be 
desirable  to  define  the  evaluation  points  while  creating  a 
primitive. 

There  are  two  successor  operators  which  generate  the 
graph  nodes  of  a  search  process.  A  start  operator  generates 
the  start  nodes  of  this  process  from  the  initial  orientation 
of  £  in  0(x,y,z,w).  These  nodes  are  put  on  OPEN  in  step  (1) 
of  algorithm  (4.21).  An  expansion  operator  generates  the 
successor  nodes  of  the  node  with  the  currently  best  spatial 
registration  in  step  (5)  of  this  algorithm. 

The  start  operator  can  use  one  of  two  approaches  in 
generating  the  initial  estimates  of  the  surface  match:  (1) 
align  known  surface  points,  or  (2)  align  orientation  of  sur¬ 
face-normal  vectors.  In  the  first  method,  there  are  L  (2  4) 
points  located  on  each  surface.  If  they  are  not  individu¬ 
ally  matched,  all  of  their  permutations  must  be  matched, 
generating  LI  initial  nodes.  This  number,  however,  can  be 


105 


substancially  reduced  if  there  is  (a)  k  priori  knowledge  of 
their  matches,  (b)  the  matches  can  be  generated  with  a 
method  such  as  stochastic  labeling  [111.  Once  the  matches 
of  individual  points  have  been  established,  the  transforma¬ 
tion  T  is  computed  from  equations  (4.6-11)  and  procedure 
(4.12)  [Section  4.21.  This  initial  estimate  of  the  spatial 
registration  T,  however,  contains  errors  caused  by  measure¬ 
ments  of  the  matched  points,  uncertainties  of  their  matches, 
and  the  orthogonality  requirement  of  the  rotation  matrix  R. 
This  registration  is,  therefore,  improved  by  a  further 
search  for  a  better  alignment. 

The  second  method  generates  the  initial  estimates  of 

the  spatial  registration  by  alignment  of  surface-normal 

« 

vectors  present  in  the  bounding  volumes  of  the  surface  quad¬ 
trees  and  representing  orientation  of  the  surface  elements 
within  the  volume  [Section  3.3).  The  two  quadtrees,  G  and 
fi',  are  traversed  from  their  root  nodes  to  bounding-volume 
nodes  of  selected  size  (volume),  and  the  normal  vectors  in 
these  volumes  are  spatially  aligned.  Since  the  four  succes¬ 
sor  nodes  in  a  quadtree  are  spatially  ordered  [Figure  3.4), 
the  two  strategies  of  matching  completely  or  partially  over¬ 
lapping  surfaces  [Section  4.2)  can  be  implemented.  If  there 
is  k  priori  knowledge  that  the  surfaces  partially  overlap, 
then  each  quadtree  is  traversed  only  to  the  bounding-volume 
nodes  near  the  parametric  boundaries  of  the  surface  descrip¬ 
tion.  However,  if  the  surfaces  completely  overlap  then  each 
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quadtree  is  traversed  to  all  the  bounding-volume  nodes  of 


the  selected  size.  Each  visited  bounding-volume  node  in  G 
is  matched  with  each  visited  bounding-volume  node  in  £'. 

The  surface-normal  vectors  in  two  bounding  volumes,  £ 
and  £' ,  of  £  and  £',  respectively,  are  aligned  as  follows: 
the  translation,  "t,  is  computed  by  translating  the  center  of 
fi  to  the  center  of  £' : 


(4.24) 


and  the  rotation,  R,  of  £  around  this  point  is  computed  by 
rotating  N  to  N 1 ,  and  two  additional  orthogonal  vectors  N-l 
and  N2  in  £  to  Nj  and  in  £' ,  so  that: 

N'  =  R  N 

N;[  =  R  Nx  (4.25) 

H  -  R  n2 

where 

N  .Nx  ■  N  *N2  *  Nx *N2  *  0 
N '  -Nj_  =  N '  •  *  0. 


These  three  pairs  of  orthogonal  unit  vectors  are  required  to 
compute  the  nine  coefficients  of  R.  After  N  has  been 
rotated  to  N',  surface  £  can  still  be  rotated  around  N  [Fig¬ 
ure  4.41.  This  rotation  is  determined  by  projecting  the 
centers  of  bounding  volumes  which  are  the  successors  of  B 
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and  B',  respectively,  to  planes  perpendicular  to  N  and  N', 
respectively.  £  is  then  rotated  around  N  to  match  every 
projected  successor  of  £  with  every  projected  successor  of 
£'  and  to  generate  up  to  four  initial  surface  registrations. 
Each  rotation  is  found  by  minimizing  the  sum  of  distances 
between  the  projected  successors  of  £  and  the  projected  suc¬ 
cessors  of  £' . 

This  second  type  the  start  operator  is  implemented  by 
the  matching  algorithm  because  it  uses  only  the  surface  in¬ 
formation  provided  by  the  quadtree  surface  representation 
[Chapter  31.  It  does  not  require  detection  of  special 
points  or  features  on  the  two  surfaces  and  their  matching  or 
any  other  additional  extraneous  processing.  If  the  size  of 
the  bounding  volumes,  aligned  by  the  start  operator,  is 
small  then  a  good  estimate  of  the  spatial  registration,  re¬ 
sulting  in  a  short  search,  is  found.  This,  however,  usually 
causes  a  large  number  of  initial  nodes  to  be  generated. 

The  expansion  operator  attempts  to  improve  the  spatial 
registration  of  the  best  node  present  on  OPEN  in  step  (5)  of 
the  matching  algorithm  (4.21)  by  generating  its  successors 
with  new  3D  rigid  transformations.  It  has  these  two 
strategies  available: 

(1)  compute  T  from  the  coordinates  of  Ti(x,y,z,w)^  and 
h'(x,y,z,w)^  for  all  jj  valid  intersection  points; 
and 

(2)  modify  values  of  the  individual  transformation  pa- 


rameters  of  T. 


The  first  strategy  generates  a  successor  whose  trans¬ 
formation  T  is  computed  by  minimizing  the  sum  of  distances 
of  all  U  points  Ti(x,y,z,w)i  and  their  intersections 
"E *  (x, y, z,w) ^  as  described  in  Section  4.2  by  equations 
(4.6-11)  and  procedure  (4.12).  Additional  successor  nodes 
may  also  be  generated  with  transformations  which  are  either 
only  translations  or  rotations.  A  translation  transforma¬ 
tion,  t,  is  computed  by  minimizing  the  sum  of  distances: 

N 

min  h'i  -  +  t  |2],  (4.26) 

i^l 

a  rotation  transformation,  R,  is  computed  by  minimizing  the 
sum  of  distances: 

t 

N 

min  ^  |  h^  -  R  |2  +  |  Rfc  | 2  ^  ,  (4.27) 

1=1 

or  by  minimizing  the  sum  of  angular  differences  of  sur¬ 
face-normal  vectors: 

N 

min  CD  R  •  N'i  -  1.0  |.  (4.28) 

i=l 

Equations  (4.26)  and  (4.27)  are  evaluated  for  the  coordi¬ 
nates  of  points  H(x,yfzfw)i  and  h'  (xf  y,  z,w)  i ;  equation 
(4.28)  is  evaluated  for  surface-normal  vectors  and  N ' ^  at 
Ti(x,y,z,w)  ^  and  h'  (x,y,z,w)ir  respectively.  Only  the  val¬ 
id  intersection  points  are,  of  course,  used  to  evaluate 
these  equations.  The  rotation  matrices  computed  from  (4.27) 


and  (4.28)  are  made  orthogonal  by  procedure  (4.12).  This 
strategy,  therefore,  may  produce  four  successor  nodes  by 
minimizing  spatial  differences  between  the  evaluation  points 
of  £  and  their  intersections  on  £'. 

The  second  expansion  strategy  modifies  the  values  of 
the  individual  transformation  parameters  of  T.  The  current 
value  of  each  active  parameter  is  incremented  and  decrement¬ 
ed  by  a  value  proportional  to  the  current  surface  difference 
D (T) .  There  are,  usually,  six  active  parameters  -  tx,  ty, 
tz,  rx,  ry,  rz  -  whose  values  are  modified,  thus  producing 
12  successor  nodes  n^  to  n^2  [Figure  4.5].  There  are  two 
functions,  At(e(n))  and  Ar(e(n)),  which  compute  modifica¬ 
tions  of  the  translation  and  rotation  parameters,  respec¬ 
tively: 

Atx(n^  )  =  +At(e<n))  |Atx(n)| 

Atx(n2  )  =  -At(e(n))  |Atx(n)| 

_  (4.29) 

•  •  •  • 

Arz(n^)  =+Ar(e(n))  |Arz(n)| 

Arz(n^2)  *  -Ar(e(n))  |Arz(n)| 

The  refinement  functions  At(e(n))  and  Arte  (n))  are  func¬ 
tions  of  the  value  of  the  evaluation  function  e  (n)  of  the 
current  node  n.  If  they  halve  the  parameter  increments  of 
the  current  node  then  this  strategy  approximates  a  binary 
search.  The  values  of  the  parameter  increments  in  a  start 
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node  are  set  to  allow  the  center  of  a  bounding  volume  B  to 
be  translated  anywhere  within  B',  and  to  allow  the  sur¬ 
face-normal  vector  ii  in  S  to  be  rotated  by  90  degrees  from 
N'  in  £' . 

A  combination  of  both  strategies  is  currently  employed 
by  the  matching  algorithm.  The  second  strategy  is  more 
effective  near  the  beginning  of  a  search;  as  the  regis¬ 
tration  improves,  the  nodes  generated  by  the  first  strategy 
converge  to  the  solution  substantially  faster  since  each 
contains  modifications  of  all  active  parameters  rather  than 
only  one  parameter.  Nodes  whose  transformations  are 
computed  only  as  translations  (4.26)  or  rotations  (4.27-28) 
are  usually  not  as  effective  as  the  nodes  whose  transforma¬ 
tions  contain  both  translations  and  rotations. 

The  matching  algorithm  has  been  able  to  find  the  cor¬ 
rect  spatial  registration  of  complicated  multiply-curved 
surface  segments,  described  by  continuous,  differentiable 
surface  representations.  If  the  registration  of  the  sur¬ 
faces  is  ambiguous  or  the  common  (overlapping)  surface  area 
is  small,  the  algorithm  finds  a  solution  and  can  be  restart¬ 
ed  to  find  further  solutions  from  nodes  still  present  on 
OPEN. 


4.4  Merging  Algorithms 


The  second  major  problem  in  assembling  complete  models 
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of  3D  objects  from  surface  segments  is  that  of  merging  two 
or  more  surface  descriptions  of  overlapping  segments,  de¬ 
fined  in  the  same  object  coordinate  system,  into  a  single 
surface  description,  or  converting  them  into  a  volumetric 
representation.  This  section  presents  several  algorithms 
which: 

(a)  merge  overlapping  surface  descriptions,  given  in 
the  same  object  coordinate  system,  into  a  single 
surface  description; 

(b)  merge  overlapping  surface  descriptions  of  a  closed 
object,  given  in  the  same  coordinate  system,  into  a 
single  volumetric  representation; 

(c)  reparametrize  a  parametric  surface  representation, 
or  parametrize  an  algebraic  surface  representation; 
and 

(d)  convert  a  surface  representation  of  a  closed  object 
into  a  volumetric  representation. 

During  the  process  of  generating  a  model  of  an  object 
from  surface  segments  we  need  to  merge  two  surface  descrip¬ 
tions  Gj'  and  G^i  under  these  two  circumstances: 

(a)  G^  and  G^i,  originally  defined  in  object  coordinate 
systems  0(x,y,z,w)|(  and  0  (x,y,  z,w)  and  parameter 
coordinate  systems  P(u,v,w)j  and  P(u,v,w)ji,  re¬ 
spectively,  have  been  matched  and  transformed  into 
0  (x, y, z, w) ^ (  and  need  to  be  merged  into  a  single 
parameter  coordinate  system  P(u,v,w)  (i.e.,  during 


the  surf ace-acquisistion  step,  the  object  was  moved 


into  a  new  orientation,  which,  in  effect,  also 
moved  the  orientation  of  the  parametrizing  grid  re¬ 
gardless  of  whether  the  grid  was  actually 
physically  moved);  or 

(b)  Gk  and  Gk»,  defined  in  the  same  object  coordinate 
system  0(x,y,z,w)k  =  0(x,y,z,w)ki  but  different  pa¬ 
rameter  coordinate  systems  P(u,v,w)j  and 
P(u,v,w)ji,  respectively,  need  to  be  merged  into 
the  same  parameter  coordinate  system  P(u,v,w) 
(i.e.,  during  the  surface-acquisition  step,  the  ob¬ 
ject  remained  in  the  same  orientation  but  the  pa¬ 
rametrizing  grid  was  moved  into  a  new  orientation). 

Finally,  when  all  the  surface  segments  of  a  closed  3D 
object  have  been  matched  and  are  defined  in  the  same  object 
coordinate  system,  it  is  also  desirable  to  merge  them  into  a 
volumetric  representation.  The  last  algorithm  in  this  sec¬ 
tion  converts  a  surface  representation,  which  consists  of  o- 
verlapping  surface  segments  or  a  single  surface  description, 
into  rectangular  parallelepiped  volumes. 

4.4.1  Transformation  of  Parameters 


G 

P 


This  method  merges 
j  and  Gji,  parametrized 
(u,v,w)j  and  P(u,v,w)^ 


two  parametric  surface  descriptions, 
by  parameter  coordinate  systems 
i,  respectively,  by  a  transformation 


[Figure  4.61 : 

u'  *  u'  <u,v)  (4.30) 

v'  *  v'  (u,v) 

which  has  the  property  that  the  functions  jj*  and  have 

continuous  partial  derivatives,  and  the  transformation  can 
be  inverted  [221.  A  linear  transformation  from  parametric 
space  P(u,v,w)j  to  parametric  space  P(u,v,w)ji  is  expressed 
as: 


(4.31) 
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The  six  coefficients  of  T{pj-»Pj,}  can  be  determined  from 
three  non-colinear  points  E(u,v,w)j  in  P(u,v,w)j  and  their 
corresponding  images  Ti(u,v,w)ji  in  P(u,v,w)ji.  They  are  ob¬ 
tained  by  rewriting  (4.31)  as  two  simultaneous  linear 
equations: 


(4.32) 
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Figure  4.6 


Surface  merging  of  two  parametric  surfaces  by 
transformation  of  a  parametric  coordinate  sys¬ 
tem:  (a)  surfaces  Gj  and  Gj ?  in  P(u,v,w)j  and 
P(u,v,w)jt,  respectively r  (d)  P(u,v,w)j  trans¬ 
formed  to  P(u,vfw)^i  merging  G^  and  G^i 


(4.33) 


15  T  =  D' 


for  each  point  E(u,v,w)j  and  its  image  h(u,v,w)ji.  For  L  2 
3  points,  system  (4.33)  is  again  solved  by  the  least-squares 
method: 


(4.34) 


The  corresponding  parametric  points  in  Gj  and  G j .  are 
again  found  by  the  ray-tracing  method.  At  each  control 
point  Ti(m,n)ji  of  G j  i ,  where  0  1  m  1  M  and  0  1  n  1  N,  a  ray 
is  cast  in  the  direction  of  the  surface  normal  vector  and 
intersected  with  Gj.  If  the  nearest  intersection  point, 
Ti(u,v)j,  is  such  that 

|  h(x,y,z,w)j  -  h(x,y,z,w)ji  |  <  epsilon  (4.35) 


where  h(x,y,z,w)j  =  h(u,v)j  and  h(x,y,z,w)ji  =  h(m,n)ji, 
then  the  parameter  values  (m,n)  and  (u,v)  belong  to  the  same 
surface  point,  h(x,y,z,w)j  =  h  (x, y, z, w) j i ,  and  form  equation 
(4.32).  If  three  or  more  such  pairs  of  points  are  found, 
then  the  transformation  T  is  computed  with  equation  (4.34). 
The  two  surfaces  [Figure  4.6(a)]  are  merged  by  expanding  G j • 
to  include  Gj.  A  new  set  of  control  points  in  the 
P(u,v,w)ji  parameter  coordinate  system  is  computed  on  a 
larger  grid  Mmin  1ml  M ^ax  and  Nmin  1  n  1  Nmax  [Figure 
4.6(b)].  A  control  point  Ti(m,n),  defined  in  G  j  • »  is 
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directly  copied  to  the  new  grid.  A  control  point  "E(m,n), 
not  defined  in  Gji,  is  transformed  by  T(Pji— »Pj)  to 
P<u,v,w)j  and  obtained  from  Gj.  If  the  point  is  also  not 
defined  in  Gjf  then  it  remains  undefined  in  the  new  grid. 
This  procedure  is  outlined  as  follows: 

(4.36) 


procedure  merge_tr ansf orm; 
begin 

[compute  transformation) 

L  :=  0; 

for  n  :=  0  step  1  until  M  do 
begin 

for  m  :=  0  step  1  until  N  do 
begin 

make_ray (ray (m, n) ,G^ i ) ; 
intersect (ray (m, n) ,  Gj , point (u, v) ) ; 
if  point (u,v)  valid  then 

L  :=  L  +  1? 
match  (m,  n,  u,  v,L) ; 
end; 

end; 

und; 

if.  L  <  3  then  return  (fail ) : 
solve  (match,  L,  T) ; 

[merge  surfaces  to  a  new  network) 

n  :=  Nmin  step  l  until  Nmax  &Q 
begin 

“min  JSiUfi  1  until  m^x  & 

begin 

find  (point  (m,n)  ,G^  • ) ; 
if  point  (m,n)  void  then 

transform (u,v,m,n,T); 
find  (point (u, v) ,Gj ) ; 
if  point(u,v)  valid  then 
point(m,n)  :=  point(u,v); 

end; 

output  (point  (m,  n) ) ; 
end; 
end; 
return: 
end; 


In  this  method  the  parameter  coordinate  system  of  G^ i 


remains  in  place  while  the  parameter  coordinate  system  of  Gj 
is  transformed,  potentially  disrupting  the  parametrization 
of  Gj.  The  next  section  contains  several  algorithms  which 
parametrize  both  surfaces  with  a  new  parameter  coordinate 
system. 

4.4.2  Projection  of  Parameters 

This  technique  merges  two  or  more  parametric  surfaces 
into  a  single  parameter  coordinate  system  by  projecting  this 
parameter  system  on  the  surfaces  [Figure  4.7].  This 
approach  is  similar  to  that  used  in  Chapter  2  to  acquire  3D 
surface  data  by  physically  projecting  an  actual  grid.  Here 
the  projection  is  done  to  surface  models  in  a  data  base. 
The  projected  parameter  coordinate  system  is  usually  defined 
by  a  plane,  and  the  projection  is  orthographic  [Figure 
4.7(a)],  or  perspective  [Figure  4.7(b)].  Alternatively,  the 
parameter  coordinate  system  could  be  defined  by  a  spherical 
or  cylindrical  surface  surrounding  the  surfaces  to  be  param¬ 
etrized.  In  addition  to  reparametrizing  parametric  sur¬ 
faces,  this  method  can  also  be  used  to  parametrize  a  surface 
defined  only  by  an  algebraic  representation. 

This  procedure  defines  a  2D  parameter  coordinate  system 
P(u,v,w)  in  a  3D  object  coordinate  system  0(x,y,z,w).  At 
sampling  distances  Au  and  Av  in  the  plane,  rays  are  cast 
from  the  plane  to  the  surfaces;  the  first  intersection  of  a 


Surface  merging  by  reparametrization:  a  param 
eter  coordinate  system  is  projected  on  sur 
faces  with  (a)  orthographic  projection,  (b 
perspective  projection 


ray  with  a  surface  then  becomes  a  control  point  on  the  newly 
parametrized  surface.  If  overlapping  surface  segments  are 
being  merged,  then  the  nearest  intersections  of  a  ray  and 


the  different  segments,  which  are  within  a  given  tolerance, 
are  averaged  into  a  single  intersection  point.  The  inter¬ 
section  point  of  a  ray  at  h(u,v)  is  connected  to  the  inter¬ 
section  points  found  by  the  four  adjacent  rays  at  Ti(u+Au), 
!i(u-Au,v),  h(u,v+Av),  and  h(u,v-Av).  The  algorithm, 
therefore,  defines  the  usual  bivariate  network  of  3D  points, 
which  is  subsequently  converted  into  a  quadtree  of  composite 
bicubic  patches: 

(4.37) 

procedure  merge_parametrize; 
bsqjn 

make_projection  (plane) ; 
lat  V  :=  vmin  Av  until  vmax  da 

begin 

igjL  U  :=  UjHiH  Step  Au  until  umax  da 
begin 

make_ray (ray (u, v) , plane) ; 
intersect  (ray (u, v) ,G, point (u, v) ) ; 
connect  (point  (u,v) , point  (u- Au,v) ) ; 
connect  (point  (u,v) , point  (u,  v- Av) ) ; 
end: 
end; 
return; 
end; 


The  algorithm  is  illustrated  in  Figure  4.8 (a-c).  A  sphere 
[Figure  4.8(a)]  is  parametrized  by  orthographic  projection 
[Figure  4.8(b)],  and  perspective  projection  [Figure  4.8(c)] 
of  a  planar  parameter  coordinate  system.  There  arc  32  x  32 
sample  rays  cast  from  the  projection  plane  to  the  sphere. 


Figure  4.8 


Parametrization  of  a  sphere  with  a  projected 
parameter  system:  (a)  image  of  a  sphere,  (b) 
orthographic  projection,  (c)  perspective  pro¬ 
jection,  (d)  closed  parametric  representation 
from  orthographic  projection 


Algorithm  (4.37)  processes  only  the  surface  points 
nearest  to  the  projected  plane  along  a  sample  ray.  For  a 
solid  object,  as  shown  in  Figure  4.8  (b-c),  it  generates  a 
parametric  surface  network  of  only  the  nearest  surface  sec¬ 
tions  facing  the  projection  plane.  However,  the  algorithm 
can  be  modified  to  produce  a  closed  parametric  network  of 
the  complete  object.  First,  all  intersections  of  a  ray  and 
the  surfaces  of  a  solid  object  are  computed.  There  must  be 
2  £  such  intersections  where  £  is  the  number  of  line  seg¬ 
ments  of  a  ray  inside  the  object.  Intersection  points  where 
a  ray  is  tangential  to  a  surface,  and  only  one  intersection 
point  is  found,  are  eliminated  by  testing  the  surface-normal 
vector,  N,  at  the  point  of  intersection.  If  N  is  nearly  or¬ 
thogonal,  within  a  tolerance,  to  the  ray  then  the  intersec¬ 
tion  is  deleted.  Intersections  of  overlapping  surface  seg¬ 
ments  are  again  averaged  into  a  single  intersection.  The 
second  modification  of  algorithm  (4.37)  affects  the  creation 
of  connections  of  the  intersection  points  located  by  adja¬ 
cent  rays:  if  a  point  on  an  adjacent  ray  is  closer  to  the 
current  point  than  the  next  point  on  the  current  ray  then  it 
is  connected  to  the  current  point;  otherwise,  the  next 
point  on  the  current  ray  is  connected  to  the  current  point. 
A  closed  network  of  3D  points  is  created  by  the  new  algo¬ 
rithm,  assuming  that  the  object  does  not  have  any  extrusions 
or  protrusions  thinner  than  the  sampling  distances  Au  and 


procedure  merge_close; 
begin 

make_projection (plane) ; 

v  :=  vmin  Av  vmax  ^ 

begin 

l&L  u  :=  step  Au  until  da 
begin 

make_ray (ray (u,  v) , plane) ; 

intersect (ray (u, v) fG, point (u,v,s),N(s),S2); 
delete (ray (u, v) , point (u,v,s),N(s),S2 ,0.010); 
for  s  :=  1  step  2  mfcil  S2  da 
begin 

if  1  point  (ur  v,  s)-point  (u,  v,  s+1)  I  > 

Ipoint  (u,v,s)-point  (u-Au,v,s)  I  then 
begin 

connect  (point  (u,  v,  s) ,  point  (u-  Au,  v,  s) ) ; 
connect  (point  (u,  v,  s+1) ,  point  (u-Au,v,s+l) ) ; 
end; 
else 

connect (point (u, v, s) , point (u, v, s+1) ) ; 
if  Ipoint  (u, v,s)-point  (u,  v,  s+1)  I  > 

I  point  (u,  v,s)-point  (u,  v-  Av,  s+1)  I  then 

begin 

connect  (point  (u,  v,  s) ,  point  (u,  v-  Av,  s) ) ; 
connect  (point  (u,  v,  s+1) ,  point  (u,  v-  Av,  s+1) ) ; 
end; 
else 

connect  (point  (u,  v,  s) ,  point  (u,  v,  s+1) ) ; 
end; 

and; 

.end; 

return; 

and; 


A  closed  network  of  3D  points  of  a  sphere  (Figure  4.8(a)], 
generated  by  this  algorithm,  is  shown  in  Figure  4.8(d). 
There  are  again  32  x  32  sample  rays  cast  from  the  param  ;ter 
plane. 


The  previous  algorithm  (4.38)  generated  a  closed 
network  of  3D  surface  points  of  a  solid  object.  While  the 
surface  sections  approximately  parallel  to  the  projection 


plane  were  closely  sampled  (with  distances  of  adjacent 
points  on  the  order  of  Au  or  Av),  the  surface  sections  ap¬ 
proximately  orthogonal  to  the  projected  plane  were  coarsely 
sampled.  An  improved  algorithm  samples  the  object  by  rays 
projected  from  three  mutually  orthogonal  planes  [Figure 
4.9],  For  each  projection  plane  only  the  intersection 
points,  where  the  angular  difference  between  a  ray  and  one 
of  the  components  of  the  surface-normal  vector  is  less  than 
45  degrees,  are  retained.  They  are  connected  to  the 
corresponding  points  found  by  the  adjacent  rays,  as  in  algo¬ 
rithm  (4.37),  and  form  disjoint  networks.  After  these 
networks  have  been  obtained  for  all  three  planes,  they  are 
connected  into  a  single  network  by  joining  points  in  differ¬ 
ent  networks  within  the  sampling  distance: 

(4.39) 

procedure  merge_orthogonal ; 
begin 

for  p  :*=  1  step  1  until  3  do 

begin 

make_projection  (plane  (p) ) ; 

£SL  V  :=  vmin(p)  S-t££  Av  until  vmax  (p) 
begin 

12JL  u  :=  u^fp)  step  Au  until  u„,ax(p)  Oa 
begin 

make_ray (ray (u, v) , plane (p) ) ; 
intersect (ray (u, v) ,G, point (u, v, s) ,N (s) ,S2) ; 
delete (ray (u, v) , point (u, v, s) ,N (s) ,S2 , 0. 707) ; 
for  s  : =  1  step  2  until  S2  do 
begin 

connect  (point  (u,  v,  s) ,  point  (u-  Au,  v,  s) )  ; 
connect  (point  (u,  v,  s) ,  point  (u,  v-  Av,  s) ) ; 
connect (point (u,  v,  s+1) , 

point  (u-  Au,  v,  s+1) ) ; 
connect  (point  (u,  v,  s+1) , 

point  (u,v- Av,s+1) ) ; 


* 


Figure  4.9  Surface  merging  by  reparametr ization:  three 

mutually  orthogonal  parameter  coordinate  sys¬ 
tems  are  or thographically  projected  on  the 
surfaces  of  a  closed  object 
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create_networks (point) ; 

return; 

end; 

This  approach  is  illustrated  with  a  sphere  [Figure  4.8(a)] 
in  Figure  4.10.  The  three  orthogonal  planes  are  projected 
to  generate  the  top  and  bottom  [Figure  4.10  (a) ],  front  and 
back  [Figure  4.10(b)],  and  left  and  right  [Figure  4.10(c)] 
networks  which  are  then  connected  into  a  single  network 
[Figure  4.10(d)].  There  were  24  x  24  sample  rays  cast  from 
each  projection  plane. 

The  last  two  algorithms  (4.38  and  4.39)  can  produce 
degenerate  surface  regions  bounded  by  three  or  six  parame¬ 
tric  curves.  Such  regions  need  to  be  specially  processed 
when  the  parametrized  surface  is  being  converted  into 
topologically  rectangular  patches.  A  region  bounded  by 
three  curves  is  handled  as  a  degenerate,  topologically 
triangular  patch.  A  region  bounded  by  six  curves  is  split 
into  two  topologically  rectangular  patches. 

In  order  that  algorithms  (4.37-39)  generate  3D  surface 
data  consistent  with  the  data  provided  by  the  photo- 
grammetric  reconstruction  method  [Section  2.4],  the  sur¬ 
face-normal  vector  is  computed  at  each  control  point  of  a 
parametric  network.  Normal  vectors  of  the  surface  of  a  sol¬ 
id  object  are  oriented  to  the  outside  of  the  object,  normal 


The  last  merging  algorithm  converts  either  a  set  of  o 


verlapping  surface  segments  or  a  single  surface  description, 
both  of  which  enclose  a  solid  object,  into  a  volume  repre¬ 
sentation  made  of  parallel  parallelepipeds.  The  algorithm 
partitions  the  object  into  parallelepipeds  by  intersecting 
parallel  rays  cast  from  a  projection  plane  with  the  surfaces 
of  the  object.  All  intersections  of  a  ray  with  the  object 
are  computed  as  in  algorithm  (4.38).  Each  line  segment  of  a 
ray  inside  the  object  defines  a  parallelepiped  centered  a- 
round  the  segment  with  a  Au  x  Av  cross-section  (Figure 
4.111 : 

(4.40) 

procedure  merge_solid; 

begin 

make_projection (plane) ; 

LSU.  v  :=  vmin+Av/2  step  Av  until  vmax  da 
begin 

lai  u  :=  Uj^+Au/2  step  Au  until  u,^  dfl. 
begin 

make_ray (ray (u, v) , plane) ; 

intersect (ray (u, v) , G, point  (u, v, s) ,N (s) ,S2) ; 
for  s  :=  1  step  2  until  S2  do 
begin 

parallelepiped  (point  (u,  v,  s) ,  point  (u,v,s+l), 
ur  Au/V,  Av,N(s),N(s+l)) ; 


return; 

end; 

The  algorithm  evaluates  the  surface-normal  vector  at  each 
intersection  of  a  ray  with  the  object  and  attaches  it  to  the 
parallelepiped.  The  normal  vector  is  oriented  to  the  out¬ 
side  of  the  solid.  This  representation  of  a  solid  object  is 
further  converted  into  the  octree  representation  144]  to 
allow  efficient  analysis  and  shaded  display. 

This  algorithm  is  illustrated  with  a  sphere  in  Figure 
4.12.  The  line  segments  of  16  x  16  sample  rays  inside  the 
sphere  are  shown  in  Figure  4.12(a).  The  16  x  16  completed 
parallelepipeds  are  shown  in  Figure  4.12(b).  These  paral¬ 
lelepipeds,  converted  to  an  octree  at  16  x  16  x  16  cube 
resolution,  are  shown  as  a  shaded  image  in  Figure  4.12(c). 
The  shading  is  computed  with  surface-normal  vectors  gener¬ 
ated  by  the  conversion  algorithm  (4.40)  from  the  algebraical 
representation  of  the  sphere  and  attached  to  each  cube.  A 
higher-resolution  octree  representation  of  the  sphere,  at  64 
x  64  x  64  cubes,  is  shown  in  Figure  4.12(d). 

The  advantages  of  using  an  analytical  model  with  a  pa¬ 
rametric  or  algebraic  surface  representation  and  converting 
it  to  a  polyhedra  representation  are:  (a)  the  analytical 
model  is  more  compact,  (b)  arbitrary  precision  of  the 
conversion  is  available,  and  (c)  analysis  and  display  of  the 
polyhedra  model  are  usually  faster. 
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(c) 


(d) 


4.12  Representation  of  a  sphere  by  parallelepi¬ 
peds:  (a)  inside  ray  segments  (16  x  16  rays), 

(b)  completed  parallelepipeds  (16  x  16  paral¬ 
lelepipeds),  (c)  low-resolution  octree  (16  x 
16  x  16  cubes),  (d)  higher- resolution  octree 
(64  x  64  x  64  cubes) 
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5 


ications  of  the  Matching  Algorithm 


AEEl 

This  section  describes  applications  of  the  matching  al¬ 
gorithm  to  these  tasks:  Cl)  surface  matching  to  build  com¬ 
plete  models  of  objects  from  surface  segments,  (2)  surface 
recognition  by  matching  a  surface  description  with  a  set  of 
objects  in  a  library  of  objects,  and  (3)  surface  matching  to 
segment  an  object  into  surface  or  volume  primitives. 

4.5.1  Object  Modeling 

Using  the  surface  acquisition  method  developed  in  Chap¬ 
ter  2  and  the  surface  matching  and  merging  techniques 
developed  in  this  chapter  we  assemble  complete  models  of  the 
measured  objects.  Each  modeled  object  is  measured  in  a 
number  of  stable  positions;  in  each  stable  position  it  is 
illuminated  in  several  parameter  positions;  in  each  stable 
and  parameter  position  it  is  photographed  from  two  or  more 
camera  positions.  The  set  of  K  stable  positions  in  which  an 
object  is  photographed  is  denoted  by: 

0(x,y,z,w)^,  ..,  0(x,y,z,w)k,  ..,  C(x,y,z,w)K. 

In  a  stable  position  j£  the  object  is  photographed  in  jlk  pa¬ 
rameter  positions  denoted  by: 

P(u,v,w)^k,  ..,  P(u,v,w)j^k,  ..,  P(u,v,w)j  ^ 

t  ,  k' 

Furthermore,  in  stable  position  Jj  and  parameter  position  j 


the  object  is  photographed  from  u  camera  positions. 


denoted  by: 


Q  (r,  S/ w)  ^  j  f  k,  Q  (r ,  s,  w)  j  f  k,  ..,  Q  ( r ,  s,  w)  j  j ,  k* 

The  total  number  of  parameter  positions,  J,  and  the  total 


number  of  images,  I,  are  given  by: 

K  K  J 

J  =  Jk,  and  I  =  ^  ^  Ijf 


(4.41) 


k=l 


k=l  j=l 


In  each  parameter  position  j  there  is  a  3D  surface  segment 
reconstructed  from  2jfk  images.  All  surface  segments,  given 
in  the  same  stable  position  Jc,  are  directly  merged  into  a 
single  surface  segment.  These  merged  surface  segments  are 
then  sequentially  matched  and  merged  into  the  complete  mod¬ 
el.  The  matching  procedure  (4.21)  is  set  to  the  "partially- 
overlapping-surfaces"  mode  [Section  4.2J  to  match  surface 
segments  of  the  same  object.  Normally,  the  surface  segments 
being  matched  are  assumed  to  have  been  obtained  from  differ¬ 
ent  stable  positions  of  the  object,  and  the  algorithm  uses 
six  active  transformation  parameters  -  tx,  hy,  tz,  rx,  rv, 
rz  -  of  the  general  transformation  (4.1).  However,  if  there 
is  4  priori  knowledge  that  the  segments  were  obtained  from 
the  same  stable  position  of  the  object  then  the  algorithm 
uses  only  three  active  parameters  -  tx,  t%,  rz  -  of  the 
limited  transformation  (4.2).  When  all  surface  segments  of 
a  solid  object  have  been  matched  they  are  merged  into  a 
single  surface  description  with  algorithm  (4.39)  and  a  vol¬ 
ume  description  with  algorithm  (4.40).  The  entire  modeling 


process  is  outlined  in  the  following  algorithm: 


(4.42) 


procedure  model ; 

.begin 

for  k  :=  1  step  1  until  K  £& 
begin 

JiSZL  j  *■  1  step  1  until  Jk  £& 
begin 

for  i  :=  1  step  1  until  I-;  u  do 
begin 

obtain  surface  information  from  image  Qi  j.  u ; 

fiM;  '3' 

reconstruct  surface  segment  j; 
if  j  >  1  then 

merge  surface  segments  j  and  j-1  into  j; 
end; 

if  k  >  1  then 
begin 

match  surface  segments  k-1  and  j; 
merge  surface  segments  k-1  and  j  into  k; 
end; 
else 
begin 

assign  surface  segment  j  to  k; 
end; 

£nd; 

merge_orthogonal ; 

merge_solid; 

return,- 

end; 


Experimental  results  of  this  modeling  process  are  illus¬ 
trated  in  Chapter  6. 


4.5.2  surface  BfigaanitiaB 


The  next  proposed  application  of  the  matching  algorithm 
is  to  surface  recognition.  Here,  the  algorithm  matches  a  3D 
surface  segment  of  an  object  obtained  from  a  single  vantage 
point  with  a  set  of  complete  models  stored  in  a  data  base. 
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The  purpose  of  such  matching  is  (1)  to  identify  the  object, 
and  (2)  to  determine  its  orientation.  The  data  base,  where 
models  are  stored,  is  organized  as  a  single  surface  graph 
[Chapter  33,  which  is  partitioned  into  groups  of  similar  ob¬ 
jects  using  the  relational-connection  nodes  [Section  3.23. 
Once  the  recognizing  procedure  determines  that  a  partition 
of  objects  may  contain  the  surface  segment,  it  visits  all 
the  objects  within  the  partition  and  attempts  to  match  them 
with  the  segment.  The  matching  algorithm  operates  in  the 
"completely-overlapping-surf aces"  mode  [Section  4.23  and 
uses  the  general  3D  transformation  (4.1) : 

(4.43) 

procedure  recognize  (SEGMENT, LIBRARY) ; 

begin 

visit  the  next  partition  R-node  in  LIBRARY; 
if  not  found  then  return  (fail ) : 
if  PARTITION  similar  to  SEGMENT  then 
begin 

visit  each  OBJECT  in  PARTITION; 

if  match (OBJECT, SEGMENT, T)  then  return (T) ? 

££d; 

recognize (SEGMENT, R-node) ; 

return  (fail) ; 

end; 


This  approach  generates  the  orientation  of  the  complete  ob¬ 
ject  from  the  surface  segment,  and  is,  therefore,  pertinent 
to  the  3D  manipulation  of  the  actual  object. 


4.5.3  Surface  and  Volume  Segmentation 


The  last  proposed  application  of  the  matching  algorithm 
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is  to  object  segmentation.  It  is  often  necessary  to  segment 
a  3D  numerical  surface  model  into  sections  that  have  common 
shape  features  or  other  geometric  properties  before  high- 
level  tasks  can  use  such  a  model.  It  is  also  desirable  to 
find  structure  relationships  among  these  segmented  sections. 
There  are  two  approaches  to  the  segmentation  process 
described  here. 

In  the  first  approach,  the  3D  models  are  segmented  into 
surface  primitives.  A  surface  primitive  is  designed  from 
one  or  more  surface  elements  of  any  available  surface  repre¬ 
sentation.  It  should  contain  meaningful  shape  properties  of 
the  model  that  is  to  be  segmented.  The  primitives  are 
designed  at  a  standard  orientation  and  scale.  The  evalu¬ 
ation  points,  where  the  shape  difference  (4.16)  between  the 
primitive  and  a  model  is  to  be  computed,  can  be  selected  at 
this  stage.  The  matching  algorithm  then  matches  the 
designed  primitives  with  the  model.  The  algorithm  uses  the 
"completely-overlapping-surf aces"  mode  and  the  general  3D 
transformation.  Since  the  primitives  are  specified  at  a 
standard  scale,  the  matching  algorithm  must  also  use  scale 
factors  as  active  transformation  parameters,  i.e.,  the 
transformation  is  not  truly  rigid  anymore.  The  primitives 
are  processed  according  to  their  priorities;  when  a  match 
of  a  primitive  and  a  model  is  found  the  surface  of  the 
primitive  is  deleted  from  the  model  and  the  process  is 
repeated  with  the  same  primitive  until  a  match  cannot  be 


found,  then  the  next  primitive  is  processed: 


(4.44) 

procedure  segment_surf ace  (SURFACE. PRIMITIVE, pmax, OBJECT) ; 

begin 

for  P  1  step  1  until  Pmax 

begin 

while  match (OBJECT, SURFACE.PRIMITIVE  (p) ,T)  ^2 

bfiflin 

delete  (OBJECT, SURFACE.PRIMITIVE (p) ,T) ; 
output (SURFACE.PRIMITIVE  (p) , T) ; 
end; 
end; 
return; 
end; 


A  union  of  the  segmented  primitives  constitutes  the 
surface  model  of  the  object. 

In  the  second  approach,  surface  models  of  solid  objects 
are  segmented  into  volumetric  primitives.  A  volumetric 
primitive  is  designed  as  a  union  of  one  or  more  surface  ele¬ 
ments  which  completely  enclose  a  3D  space.  The  matching  al¬ 
gorithm  spatially  registers  the  surfaces  of  the  primitive 
with  the  surfaces  of  the  object.  At  each  evaluation  point 
of  a  primitive,  the  half-space  (inside  or  outside),  where 
surface  difference  (4.16)  is  to  be  computed,  is  specified. 
Therefore,  a  volumetric  primitive  can  be  spatially  register¬ 
ed  to  (a)  enclose  completely,  (b)  be  enclosed  completely, 
(c)  be  completely  outside,  or  (d)  be  partially  outside  the 
solid  object.  The  union,  intersection,  and  difference  set 
operators  can  be  used  to  construct  a  CSG  tree  [Section  3.2) 
from  the  segmented  primitives.  The  matching  algorithm  uses 
the  "partially-overlapping-surf aces"  mode  and  the  general  3D 
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transformation  with  scale  factors.  The  volume- segmentation 
process  is  similar  to  algorithm  (4.44): 


(4.45) 

procedure  segment_volume  (VOLUME_ PRIMITIVE, P  _x, 

SOLID_OBJECT) ; 

begin 

CSG_TREE  : =  empty; 

P  :=  1  .S-fc.e.P  1  until  Pmax  dfi 
begin 

:<hile  match  (SOL ID_OBJECT,  VOLUME. PRIMITIVE  (p)  ,  T)  <& 

fefiflln 

delete (SOL ID_OBJECT, VOLUME. PRIMITIVE (p) , T) ; 
interference (VOLUME_PRIMITIVE (p) ,T,SOLID_OBJECT, 
operator); 

attach (VOLUME_ PRIMITIVE (p) , T, operator , CSG_TREE) ; 
end: 
end: 

return (CSG_TREE) ; 

snsLi 


The  union-set  operator  is  applied  to  primitives  com¬ 
pletely  inside  the  solid;  the  intersection  operator  is 
applied  to  primitives  partially  inside  the  solid;  and  the 
difference  operator  is  applied  to  primitives  completely 
outside  the  solid. 


4.6  Ulus.ti.fl.fcig ns 


The  first  example  illustrates  the  matching  algorithm  by 
matching  two  surface  segments  [Figure  4.131.  Surface  seg¬ 
ment  G^  is  defined  in  object  coordinate  system  0(x,y,z,w)i 
[Figure  4.13(a)],  and  surface  segment  G2  is  defined  in  ob¬ 
ject  coordinate  system  0(x,y,z,w)2  [Figure  4.13(b)).  Sur¬ 
face  segment  G^  has  been  spatially  registered  with  surface 


Figure  4.13 


Surface  matching:  (a)  surface  segment 
0  (xf y, z, w) (b)  surface  segment  G 
0  (xf  y,  z, w)  ^  i ,  (c)  matched  surface  segme 


nu  vj i 

segments 


segment  G2  and  transformed  to  0(x#y,z,w)2  [Figure  4.13(c)]. 
The  rigid  3D  transformation  which  aligned  the  two  segments 
is: 


0.235 

0.^99 

-0.323 

15.817 

0.723 

0.213 

0.658 

15.043 

0.650 

-0.532 

-0.542 

69.560 

0.000 

0.000 

0.000 

1.000 

(4.46) 


The  matching  algorithm  generated  168  start  nodes  by  aligning 
surface  normal  vectors  in  9  bounding  volumes  of  with  sur¬ 
face  normal  vectors  in  6  bounding  volumes  of  G2.  These 
start  nodes  were  expanded  into  449  nodes  before  the  solution 
was  found. 

The  next  example  illustrates  merging  of  surface  seg¬ 
ments  using  the  second  merging  algorithm  which  projects  a 
new  parameter  coordinate  system  [Figure  4.14].  There  are 
three  surface  segmentsf  G^,  G2*  and  G3f  defined  in  the  same 
object  coordinate  system  0(x,y,z,w)lr  and  three  different 
parameter  coordinate  systems  P(ufv,w)^i'  P(u,v,w)2fif  and 
P(ufv,w)3^  [Figure  4.14(a)].  A  square  (Au  =  Av)  parame¬ 
ter  coordinate  system  was  projected  on  these  segments  with 
orthographic  projection  resulting  in  a  single  surface  seg¬ 
ment  [Figure  4.14(b)]. 

The  last  example  illustrates  the  merging  algorithm 
(4.40)  which  converts  a  surface  representation  of  a  solid 
object  into  a  volume  representation  [Figure  4.15].  The  sur- 
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face  model  of  the  J79  turbine  blade  [Figures  3.10,  3.11  and 
3.153  was  partitioned  into  parallelepiped  volumes  by  inter¬ 


secting  the  model  with  orthogonally  projected  rays.  A 
sampling  resolution  of  40  x  160  rays  converted  the  blade 
into  2706  parallelepipeds.  The  line  segments  of  each  ray 
inside  the  blade  are  shown  in  Figure  4.15(a).  The  parallel¬ 
epipeds  were  further  converted  into  the  octree  representa¬ 
tion.  The  octree  model  shown  in  Figure  4.15(b)  contains 
35808  nodes,  of  which  6132  nodes  are  full.  The  shading  was 
computed  using  the  surface-normal  vectors  obtained  from  the 
analytical  model. 

4.7  Summary 

A  surface  matching  algorithm  which  spatially  registers 
3D  surfaces  was  described  in  this  chapter.  The  algorithm, 
in  conjuction  with  several  merging  algorithms,  is  used  to 
assemble  complete  models  of  3D  objects  from  their  multiple 
surface  segments.  This  process  is  illustrated  in  Chapter  6. 
Suggested  modification  of  this  algorithm  for  object 
recognition,  or  surface  and  volume  segmentation  were  also 
presented. 


CHAPTER  5 


EXPERIMENTAL  METHODS 

This  chapter  describes  the  hardware  used  and  the 
software  developed,  in  the  course  of  this  work,  to  generate 
the  illustrative  examples  shown  in  the  three  previous  chap¬ 
ters  and  the  complete  models  that  will  appear  in  the  next 
chapter. 

5.1  Hardware 


The  photogrammetric  process,  described  in  this  report, 
required  a  camera  calibration  stand,  a  projector,  a  camera 
and  a  film  scanner.  A  camera  calibration  stand  w^s  built 
from  three  sheets  of  glass  placed  perpendicularly  to  each 
other.  The  glass  sheets  defined  the  object  coordinate  sys¬ 
tem  0(x,y,z,w),  in  which  all  surfaces  were  measured.  There 
were  ten  30  x  30  mm  calibration  marks  placed  in  the  stand, 
five  each  in  the  x  »  0  and  y  =  0  planes,  all  clearly  visible 
to  the  camera.  The  corners  of  these  marks  were  detected  in 
the  generated  images  and  used  to  compute  the  camera  trans¬ 
formation  matrices  T{0— »Q).  The  size  of  the  stand  and  the 
placement  of  the  marks  allowed  measurement  of  objects  up  to 
300  x  300  x  300  mm  in  size.  A  35  mm  slide  projector  illumi¬ 
nated  the  camera  calibration  stand  with  a  pattern  of  isopa- 
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rametric  lines  defined  in  a  parameter  coordinate  system 
P(u,v,w).  The  pattern  was  recorded  by  a  photo-densitometer 
on  a  35  mm  slide.  There  are  five  20-micron  wide  lines  per 
mm  in  both  jj  and  ^  directions  of  the  pattern.  These  lines 
cover  only  a  20  x  20  mm  area  of  the  24  x  36  mm  slide  to 
allow  normal  illumination  of  the  camera  calibration  marks. 
A  35  mm  camera  with  a  55  mm  focal-length  lens  was  used  to 
photograph  the  measured  surfaces  on  a  high-contrast,  black- 
and-white  film  with  high  resolving  power.  These  images  were 
digitized,  in  Q(r,s,w)  image  coordinate  systems,  by  the 
photo-densitometer  which  sampled  the  film  at  20-micron 
intervals  (50  samples  per  mm)  in  both  £  and  s.  directions.  A 
16-bit  film  density  value  was  obtained  at  each  sampled 
point.  In  these  images,  all  the  relevant  information  (illu¬ 
minated  surfaces  and  calibration  marks)  is  contained  within 
a  24  x  30  mm  area  of  the  24  x  36  mm  image  frame.  This  means 
that  each  image  was  digitized  to  a  1200  x  1500  16-bit  pixel 
resolution. 

All  the  data  processing  was  done  by  a  32-bit 
midi-computer  with  (what  appears  to  be  at  this  time)  a 
virtually  unlimited  virtual  memory  and  a  large  disk  space. 
A  vector  drawing  and  a  color  raster  terminals  were  used  for 
display  of  graphical  and  image  data  generated  by  this  work. 
The  vector  drawing  terminal  can  display  up  to  8000  2D 
vectors  or  5500  3D  vectors  from  its  refresh  memory.  Ortho¬ 
graphic  projections  of  the  3D  vectors,  following  scale,  ro- 


tation  and  translation  transformations,  are  computed  by  the 
terminal  in  real  time,  using  special  hardware.  A  3D  surface 


description  can  thus  be  displayed  while  smoothly  "tumbling" 
on  the  screen,  thereby  increasing  the  viewer's  perception  of 
its  3D  shape.  A  hard  copy  of  the  image  displayed  on  the 
screen  can  be  made  by  an  electrostatic  printer-plotter  with 
resolution  of  about  8  binary  dots  per  mm  (see,  for  example. 
Figures  3.8,  3.10,  and  3.12).  The  raster  terminal  displays 
a  digital  color  image  from  a  512  x  512  x  24  bit  frame 
buffer.  It  was  used  here  to  display  windows  in  the 
digitized  black-and-white  images  of  the  measured  surfaces  as 
well  as  synthetic  images  computed  from  the  3D  surface  mod¬ 
els.  Although  the  synthetic  images  cannot  be  computed  in 
real  time,  the  viewer's  perception  of  the  3D  scene  shown  in 
a  single  frame  can  be  augmented  by  the  use  of  color  shading, 
shadows  and  textures.  A  black-and-white  and  a  color  hard¬ 
copy  units  are  attached  to  this  terminal  (see,  for  example. 
Figures  3.9,  3.11,  3.13). 


5.2  software 


j  In  addition  to  the  development  of  the  surface  modeling 

process,  an  important  secondary  goal  of  this  work  was  that 
the  design  and  implementation  of  the  process  be  achieved 
j  with  good-quality  software.  There  were  two  key  objectives 

in  designing  the  software: 
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(1)  flexibili tv  -  to  allow  easy  modification  and 
expansion,  since  the  system  is  primarily  experimen¬ 


tal  ;  and 

(2)  usability  -  to  allow  large  quantities  of  various 
types  of  data  to  be  processed  without  increased 
effort. 

The  software  was  implemented  as  four  command  proc¬ 
essors.  The  first  three  processors  contain  the  modeling 
procedures  described  in  Chapters  2,  3,  and  4,  respectively. 
The  fourth  processor  computes  synthetic  images  from  3D  sur¬ 
face  data  generated  by  the  modeling  process  and  from 
additional  parameters  which  enhance  the  3D  scene  (light  il¬ 
lumination,  shadows,  textures),  and  define  the  viewing 
camera  and  its  optical  system. 

A  command  processor  is  an  interactive  software  program 
which  allows  a  user  to  enter,  process,  and  store  data  with  a 
set  of  commands  and  their  parameters.  It  creates  its  own 
software  environment,  which  is  required  for  the  given 
application  and  is  separate  from  the  operating  system. 
Unlike  an  ordinary  program  which  has  its  pre-defined  input, 
processing,  and  output  parts,  a  command  processor  can 
execute  a  flexible  sequence  of  commands  appropriate  for  the 
data  being  processed.  In  the  context  of  our  experimental 
systems,  it  permits  us  to  display  the  results  of  each 
computational  step,  and  if  necessary  to  repeat  the  step  with 
different  parameters. 
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The  implemented  processors'  names  and  tasks  are  as 


follows: 

SURE  -  reconstructs  3D  surface  descriptions  from 
multiple  images  of  the  surface  illuminated  by  iso¬ 
parametric  lines; 

sued  -  converts  3D  surface  descriptions  into  hier¬ 
archical  representations,  enters  and  edits  sur¬ 
faces,  and  attaches  surface  properties; 

ALIGN  -  aligns  two  surface  representations  into  a 
common  3D  object  coordinate  space,  and  optionally 
merges  parametric  surface  representations  into  a 
common  2D  parameter  coordinate  space; 

STRAW  -  generates  shaded  synthetic  images  of  sur¬ 
face  descriptions. 

The  primary  motivations  in  the  software  design  are 
listed  below: 

(a)  Common  data  structures 

There  is  a  single  definition  of  all  the  data  struc¬ 
tures  and  parameters  within  each  processor;  this 
definition  applies  to  all  the  routines  in  the  proc¬ 
essor.  When  a  new  routine  is  added  to  a  processor 
the  data  base  definitions  are  simply  inserted  into 
the  new  module.  Because  the  processors  pass  data 
in  the  same  format  to  each  other,  there  are  large 
parts  of  the  data  structures  also  shared  among  dif¬ 
ferent  processors. 


An  abbreviation  facility  [55],  common  to  all  the 
processors,  allows  every  non-numeric  string  to  be 
abbreviated  to  the  leading  non-ambiguous  substring. 
The  abbreviation  facility  uniformly  applies  to  all 
command  names,  keywords,  directory  and  file  names, 
and  symbolic  data  identifiers. 


The  developed  software  modules  are  grouped  into 
these  categories: 

(1)  control  modules  -  control  execution  of  a  proc¬ 
essor,  contain  command  processing; 

(2)  execution  modules  -  execute  data  processing 
commands ; 

(3)  utility  modules  -  provide  various  utility  func¬ 
tions  of  a  processor; 

(4)  image  processing  modules  -  form  a  library  of 
routines  for  access  and  operations  on  2D  image 
data; 

(5)  geometrical  processing  modules  -  form  a  library 
of  routines  for  access  and  operations  on  2D  or  3D 
geometrical  data; 

(6)  input-output  modules  -  interface  to  the  data¬ 
base  management  subsystem. 

The  following  sections  describe  the  individual  command 
processors  and  any  associated  software  in  more  detail; 
Appendix  D  contains  the  syntax  of  the  actual  commands. 
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SURE  is  the  surf ace- reconst ruction  processor.  It 
implements  the  techniques  described  in  Chapter  2.  The  input 
to  this  processor  is  a  set  of  digitized  images  of  the 
measured  surfaces.  The  surfaces  in  all  such  images  are  in 
the  same  object  coordinate  system  0(x,y,z,w)^  and  are  param¬ 
etrized  by  the  same  parametric  coordinate  system  P(u,v,w)j. 
The  output  of  the  processor  is  a  list  of  reconstructed 
control  points  or  boundary  curves  of  parametric  patches. 
This  list  is  passed  to  the  surface  editor  for  conversion  of 
the  surface  information  to  surface  patches,  and  hierarchical 
structuring  to  surface  quadtrees.  The  individual  commands 
of  this  command  processor  are  listed  in  Appendix  D.l. 


SUED  is  the  .surface  suitor.  Its  main  task  is  to 
convert  a  list  of  unstructured  patch  control  points  or 
boundary  curves  into  a  structured  quadtree  of  patches, 
bounding  volumes  and  normal  vectors,  which  can  be  used  by 
the  surface-alignment  process.  In  addition,  other  surface 
representations  (planar,  quadric  as  well  as  bicubic)  can  be 
entered  by  SUED  commands  into  its  data  base.  These  repre¬ 
sentations  can  later  be  used  as  3D  shape  primitives  in  the 
alignment  process.  Various  surface  attributes  can  be 
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specified  for  each  surface  element  although  only  reflection 
and  transmission  coefficients,  and  texture  functions  are 


useful  at  the  present  time.  SUED  generates  a  structured 
surface  file  which  can  be  used  either  in  surface  alignment 
or  in  generation  of  synthetic  images.  The  commands  of  this 
processor  are  described  in  Appendix  D.2. 


5.2.3 


ALIGN  is  a  processor  which  performs  the  3D  alignment  of 
surface  descriptions,  using  the  matching  and  merging  algo¬ 
rithms  given  in  Chapter  4.  Two  surface  descriptions  are 
spatially  registered  with  algorithm  (4.21).  If  the  two  sur¬ 
faces  are  segments  of  an  object  they  are  merged  into  a 
single  parameter  coordinate  system  following  the  spatial 
registration.  After  the  two  surface  descriptions  have  been 
matched,  and  possibly  merged,  they  are  returned  back  to  SUED 
as  a  linear  list  of  surface  patches.  SUED  then  restructures 
them  into  a  single  hierarchical  description.  ALIGN  also 
performs  all  the  merging  algorithms  of  Section  4.4:  parame¬ 
tric  surfaces  are  merged  by  transformation  of  parameters; 
algebraic  and  parametric  surfaces  are  parametrized  or  repa¬ 
rametrized  by  parameter  projections;  and  surfaces  of  a  sol¬ 
id  object  are  converted  to  a  volume  representation  made  of 
parallelepipeds.  A  description  of  the  individual  commands 
of  this  processor  may  be  found  in  Appendix  D.3. 


5.2.4 


STRAW  is  a  processor  which  paints  pretty  pictures  in 
highly  non-linear  time.  It  renders  images  of  complex  3D 
scenes  with  surface-to-surface  reflections,  transparent  sur¬ 
faces  with  refractions,  diffuse  and  specular  reflections, 
image  and  texture  surface  mappings,  illumination  by  point- 
light  sources  and  shadows.  STRAW  uses  a  ray-tracing  algo¬ 
rithm  to  compute  visible  surfaces  and  a  recursive  shading 
algorithm,  developed  by  Whitted  [72],  to  compute  light 
intensity  of  the  visible  surfaces. 

A. raster  image  computed  by  this  program  is  considered 
to  be  an  array  of  square  or  rectangular  pixels.  The  3D 
sceTie  is  sampled  at  the  four  corner  points  of  each  pixel. 
In  perspective  projection  a  3D  sample  ray  from  an  image 
sample  point  to  the  center  of  projection  is  extended  into 
the  scene  and  intersected  with  the  nearest  surface  element. 
In  orthographic  projection  each  3D  sample  ray  is  perpendicu¬ 
lar  to  the  image  plane.  The  intensity  of  the  sample  point 
is  evaluated  from  visual  properties  of  the  intersected  sur¬ 
face  element  as  well  as  from  intensity  information  provided 
by  additional  rays  which  are  bounced  from  the  intersection 
point  in  the  reflection,  refraction,  and  light-source 
directions.  Information  about  each  reflected  and  refracted 
ray  is  maintained  in  a  node  of  a  binary  shading  tree. 
Finally,  the  intensity  of  a  pixel  is  computed  by  averaging 
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its  four  sample  values.  Anti-aliasing  is  performed  by 
recursively  subdividing  a  pixel,  which  has  a  large  differ¬ 
ence  in  its  sampled  values,  into  2x2  regions  and  repeating 
the  sampling  process  at  the  four  corners  of  each  new  region. 
The  intensity  of  each  region  is  again  the  average  of  its 
four  samples  and  the  final  intensity  of  the  pixel  is  the  sum 
of  all  regional  intensities,  each  weighted  by  its  area. 
This  subdivision  is  done  only  within  pixels  which  contain 
sharp  intensity  changes,  typically  caused  by  edges,  silhou¬ 
ettes  or  textures. 

STRAW  operates  on  the  structured  surface  descriptions 
provided  by  SUED.  The  ambient  intensity,  i.e.,  the  color  of 
a  surface  element  can  be  specified  as  a  constant  value  [Fig¬ 
ure  5.1(a)],  the  parameters  u,  v  of  the  element  can  be  used 
to  look  up  the  value  in  a  color  image  [Figure  5.1  (c,d)],  or 
the  value  can  be  computed  by  linear  interpolation  of  colors 
along  3D  vectors  stored  in  a  paint  table  [Figure  5.1(b)].  A 
texture  can  be  added  to  a  surface  element  with  Blinn's 
texturing  technique  [12].  This  technique  perturbs  the 
normal  vector  to  a  surface  point  in  a  direction  specified  by 
the  partial  derivatives  of  a  2D  texture  function.  Here,  the 
source  of  the  texture  function  is  a  monochrome  image  and  the 
parameters  jjr  2  of  a  surface  element  are  used  to  look  up 
values  in  the  image  from  which  the  partial  derivatives  are 
computed.  This  texturing  is  quite  effective;  however,  it 
remains  only  a  shading  illusion  -  the  silhouettes  of  sur- 
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Figure  5.1 


Examples  of  synthetic  images  generated  by 
STRAW:  (a)  Microtubular  doublet,  (b)  IPL  logo, 
(c)  Recursive  box,  (d)  Galaxy  far  far  away  ... 


faces  are  still  smooth  and  the  textures  cannot  cast  shadows. 


In  addition  to  generating  a  raster  image  with  a  pin¬ 
hole  camera  model,  STRAW  can  optionally  save  information 
about  each  image  sample  point  from  the  nodes  of  its  shading 
tree.  These  sample  points  may  be  later  converted  into  an 
actual  raster  image  by  a  post-processor.  A  post-processor, 
FOCUS,  was  designed  to  implement  the  focus  and  aperture 
functions  of  a  camera's  lens  [57,58].  It  generates  synthet¬ 
ic  images  which  are  focused  and  have  a  depth  of  field  [Fig¬ 
ures  5.2-31.  For  each  image  sample  point,  this  post-proc¬ 
essor  computes  a  point-spread  function  whose  size  and 
intensity  distribution  depend  on  the  sample's  depth  (i.e., 
the  distance  of  the  visible  surface  along  the  camera's 
optical  axis),  the  focus  distance  of  the  lens,  and  the 
aperture  size.  The  point-spread  function  is  then  convolved 
with  the  sample  point  in  the  spatial  domain.  Since  FOCUS 
has  available  information  about  each  ray  in  each  shading 
tree,  it  can  properly  defocus  reflected  [Figure  5.21  and 
refracted  [Figure  5.3]  rays  of  light.  A  post-processor, 
BLUR,  adds  motion  blur  1591,  due  to  a  finite  exposure  time 
of  real  cameras,  to  moving  surfaces  [Figure  5.4],  This  is 
accomplished  by  convolving  all  image  samples  which  belong  to 
a  moving  surface  with  a  point-spread  function  computed  from 
the  path  and  velocity  of  the  motion,  and  the  duration  of  the 
exposure.  This  convolution  can  be  performed  equally  well  in 
the  spatial  or  the  frequency  domain. 
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Opaque  carafe  (a)  made  with  pin-hole  camera 
model f  (b)  focused  on  the  carafe  and  aperture 
set  to  f/1.4 


Spheres  falling  on  a  parabolic  path:  (a)  made 
with  pin-hole  camera  model,  (b)  made  with 
motion-blur  camera  model 


Solid  objects  can  be  converted  by  STRAW  to  the  paral¬ 
lelepiped  representation  (Section  4.4.3)  with  algorithm 
(4.40).  The  line  segments  inside  a  solid  object  of  each 
ray,  cast  by  the  orthographic  projection,  are  determined  by 
the  ray-tracing  traversal  procedure.  The  surface  ambient 
intensity,  and  the  surface-normal  vector,  oriented  to  the 
outside  of  the  object,  are  evaluated  at  both  ends  of  each 
line  segment.  The  ambient  intensity  may  be  a  constant,  im¬ 
age,  or  paint- table  value;  the  surface-normal  vector  may  be 
perturbed  by  a  texture  function.  This  description  of  solid 
objects  is  passed  to  an  octree  modeling  system  (441,  which 
converts  it  to  octree  representation  for  efficient  analysis 
and  display. 

STRAW  has  also  been  successfully  used  to  generate 
stereo  pairs  of  images  that  can  be  viewed  as  transparencies 
in  a  stereo  viewer.  Command  STEREO,  given  a  vantage  point, 
a  view  point  and  a  stereo  separation  angle,  computes  the 
complementary  vantage  point.  Any  part  of  the  STRAW  data 
base  may  be  modified  between  the  generation  of  consecutive 
image  frames  that  may  be  used  in  an  animated  sequence.  This 
allows,  for  example,  movement  of  the  camera,  lights,  and 
surfaces  as  well  as  modifications  of  surface  shapes  and 
visual  properties.  The  individual  commands  of  this  software 
processor  are  listed  in  Appendix  D.4. 
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CHAPTER  6 


EXPERIMENTAL  RESULTS 

This  chapter  contains,  in  the  first  section,  an  example 
of  building  a  complete  computer  model  of  a  3D  object  by  the 
photogrammetric  reconstruction  and  the  match-and-merge 
alignment  obtained  from  several  stable  positions  of  the  ob¬ 
ject.  This  example  exercises  the  entire  modeling  process  of 
this  work.  The  second  section  shows  two  symmetric  objects 
that  were  partially  reconstructed  by  the  photogrammetric 
procedure,  and  then  completed  into  generalized  cylinder-like 
shapes  with  a  priori  knowledge  of  their  symmetry.  Each  of 
these  two  objects  is  modeled  as  a  single  sheet  of  surface 
patches.  All  the  examples  given  in  this  chapter  contain 
multi-curved  surfaces. 

6.1  Car  Model 

The  first  example  demonstrates  the  entire  modeling 
process  of  obtaining  a  3D  description  of  an  object.  It  is  a 
balsa  model  of  a  car  meant  for  the  late  1980's.  The  comput¬ 
er  model  was  computed  from  36  views  of  the  3D  object.  The 
object  was  positioned  in  six  different  object  coordinate 
systems  0(x,y,z,w)2  to  0<x,y,z,w)g.  In  each  object  coordi¬ 
nate  system  0 (x,y, z,w) it  was  consecutively  illuminated  by 
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three  parametric  coordinate  systems  P(u,v,w)^^  to  P(u,v, 
w*3/k*  While  illuminated  by  a  parametric  coordinate  system 
P (u, v,w) j ^  two  images  Q (r , s,w) ^ j ^ k  and  Q (r , s ,w) 2f j f ^  were 
recorded.  The  tree  structure  of  this  modeling  process  is 
given  in  Figure  6.1  with  i  =  36,  J.  =  18,  and  K  =  6.  Four 
low-resolution  (256  x  256  pixel)  digital  images  of  the 
original  object  are  shown  in  Figure  6.2.  A  proof  (contact) 
print  of  the  36  images  is  given  in  Figure  6.3.  Each  column 
contains  six  images  of  the  object  in  one  stable  position; 
each  pair  of  images  within  a  row  contains  the  object  param¬ 
etrized  by  the  same  parameter  system.  From  the  36  images, 
there  were  18  3D  surface  segments  reconstructed  [Figure 
6.4J .  Since  every  three  segments  were  defined  in  the  same 
object  coordinate  system,  they  were  simply  merged  into  the 
same  parametric  coordinate  system  [Figure  6.5],  In  this 
figure,  each  drawing  in  the  left  column  shows  three  surface 
segments  in  0(x,y,z,w)^  before  they  were  merged,  and  the 
corresponding  drawing  in  the  right  column  shows  the  merged 
surface  segment.  Algorithm  (4.37)  was  used  to  reparametrize 
the  surface  segments  with  orthographic  projection  and  5.0  mm 
sampling  distance.  At  last,  the  six  merged  surface  segments 
were  aligned  into  a  common  object  coordinate  system.  This 
alignment  process  consisted  of  five  iterations;  the  coordi¬ 
nate  system  of  the  first  surface  segment  remained  stationary 
while  the  other  five  segments,  one  at  a  time,  were  matched 
with  the  existing  surface  description  [Figure  6.6],  Four 


?',£■*  *14**  4>*\ 


0(xryf  ZrWjjL 


"Si**; 

SSsSiS  ,..- 
alsSffiS;:: 
SSiSfiSi?:: 


gIiSS»S55'pi 

P^isss-sr' 

,-.,.1  ,,.,J 

,|llll!>,1i!.  -5 1 

■Sl-lil 

*jS&X 


Us! 


0  (x,  y,  z,w)  2 


Figure  6.5  Surface  segments  of  Figure  6.4  merged  into  six 
segments 
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views  of  the  aligned  surface  segments  are  shown  in  this  fig¬ 
ure  together  with  the  axes  of  the  object  coordinate  system 
of  the  first  segment.  A  single  surface  model  of  the  object 
was  therefore  obtained,  and  structured  into  a  hierarchical 
description. 


This  example  illustrates  the  reconstruction  of  two  ob¬ 
jects  which  are  symmetric  around  an  axis.  Each  of  the  ob¬ 
jects  -  a  wine  carafe  and  a  coke  bottle  -  was  measured  in  a 
Single  object  coordinate  system  0(x,y ,z ,w) j,  while  param¬ 
etrized  with  a  single  PCurVfW)-^^  parametric  coordinate  sys¬ 
tem,  and  photographed  '.n  two  image  coordinate  systems: 
Q(r,s,w)^^1^2  and  Q(r,s,w) 2, 1, i*  The  section  of  each  ob¬ 
ject,  reconstructed  from  the  two  images,  was  then  completed 
into  a  generalized  cylinder-like  shape  made  of  a  single 
sheet  of  bicubic  patches.  It  was  assumed  that  each  object 
had  a  circular  cross-section  in  x-y  planes.  Intersection 
curves  of  the  reconstructed  3D  section  with  x-y  planes 
equally  spaced  between  the  ground  plane  and  the  top  of  the 
object  were  computed  and  circles  were  fitted  into  the  curves 
with  a  least-squares  error  method.  Figure  6.7  shows 
medium-resolution  (512  x  256  pixel)  digital  images  of  the 
original  objects.  Figures  6.8-6.11  show  the  results  of  this 
modeling  procedure.  Figures  6.8  and  6.10  contain  the  two 
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original  images  of  the  parametrized  objects.  Figures  6.9(a) 
and  6.11(a)  are  line  drawings  of  the  quadtree  representa¬ 
tions  of  the  two  objects  and  Figures  6.9(b)  and  6.11(b)  are 
shaded  synthetic  images  of  the  two  models.  The  point-light 
source  in  each  of  these  images  is  placed  approximately  in 
the  same  location  at  the  lamp  illuminating  the  original  ob¬ 
jects  in  Figure  6.7.  Notice  (or  rather  do  not  notice)  the 
loss  of  detail  in  the  surface  shape  visible  along  the  sil¬ 
houette  of  each  object.  This  is  partially  caused  by  the 
resolution  of  the  projected  parameteric  lines,  and  also  due 
to  using  circular  cross-sections  in  the  horizontal  direction 
and  cubic  spline  in  the  the  vertical  direction  of  the  ob¬ 
jects.  Figure  6.12  contains  two  images  of  these  models 
placed  in  a  model  of  the  camera  calibration  stand.  The  mod¬ 
el  of  the  carafe  is  also  pictured  in  Figures  5.2,  and  5.3. 


CHAPTER  7 


RETROSPECTS  AND  PROSPECTS 

This  chapter  summarizes  the  presented  work,  lists  its 
major  contributions  and  then  proposes  tasks  for  further 
research. 

7.1  Recapitulation 

The  main  objective  of  this  research  was  to  investigate 
techniques  for  computer  generation  and  processing  of  three- 
dimensional  surface  models  of  existing  physical  objects. 
The  obtained  numerical  models  were  to  be  reasonably  precise, 
complete,  and  hierarchicaly  structured.  Models  of  man-made 
objects,  applicable  to  both  computer  vision  and  CAD/CAM, 
were  of  especial  interest. 

An  image  processing  technique,  which  digitizes  3D  sur¬ 
faces  in  a  controlled  environment  (position  and  illumination 
of  the  surfaces),  was  developed  to  obtain  3D  surface  data 
located  on  a  2D  parametric  grid.  These  data,  computed 
originally  at  a  high  resolution,  are  structured  to  provide 
hierarchical  surface  descriptions.  Depending  upon  the 
specific  amount  of  surface  detail  required,  various  levels 
of  the  hierarchical  structure  are  employed.  A  matching  al¬ 
gorithm  uses  these  hierarchical  representations  to  perform 
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3D  alignment  of  3D  surfaces  which  share  common  (overlapping) 
sections. 

A  modeling  process  which  assembles  complete  models  of 
arbitrarily-shaped  3D  objects  was  developed  from  these 
capabilities.  Individual  3D  descriprions  of  the  surfaces  of 
the  object  are  first  obtained  for  different  positions  of  a 
parameter  coordinate  system  and  also  for  different  position 
of  the  object  in  its  object  coordinate  system.  The  descrip¬ 
tions,  given  in  the  same  object  system  but  different  parame¬ 
ter  systemsr  are  merged  into  a  common  parameter  system. 
Those  descriptions,  that  are  given  in  different  object  sys¬ 
tems,  are  first  spatially  matched  into  a  common  object  sys¬ 
tem,  and  then  merged  into  a  common  parameter  system.  If  3D 
data  covering  all  the  surfaces  of  an  object  are  available, 
this  process  can  iteratively  build  a  single  complete  model. 
The  model  representation  developed  here  is  versitile  enough 
to  be  directly  used  in  a  number  of  applications;  it  can 
also  be  easily  expanded  and  converted  to  other  representa¬ 
tions.  Two  immediate  applications  are  (a)  conversion  to 
octree  representation  for  fast  object  analysis  and  display 
[44],  and  (b)  generation  of  characteristic  views  for  object 
recognition  [17].  A  single  ray-tracing  procedure  provides 
all  required  information  about  a  surface  description  to  all 
application  programs.  All  geometric  operations  which  depend 
on  surface  representation  are  confined  to  this  procedure. 

The  surface  and  object  models  can  have  assigned  reflec- 


tion,  transmission,  and  texture  characteristics,  and  be 
viewed  in  computer-generated  raster  images.  The  circle  of 
computer  graphics  and  image  processing  is  then  completed  by 
applying  image  processing  techniques  to  these  synthetic  im¬ 
ages.  The  inverse  process  of  "image  restoration,"  used  in 
image  processing  to  remove  noise  from  natural  images,  is 
applied  here  to  add  noise,  caused  by  the  optical  system  of 
the  simulated  camera,  to  the  computer  graphics  images. 

7.2  Contributions 

The  following  may  be  considered  as  the  major 
accomplishments  of  this  work: 

(1)  Development  of  a  3D  surface-data  acquisition  system 
based  on  spatially-controlled  illumination  of  the 
measured  surfaces  in  multiple  photographic  images; 
reconstruction  of  surface  information  from  2D  per¬ 
spective  projections  into  a  3D  representation. 

(2)  Hierarchical  structuring  of  the  3D  surface  repre¬ 
sentation  into  surface  quadtrees  of  varying  levels 
of  detail. 

(3)  An  algorithm  for  computing  intersections  of  a  3D 
line  and  a  3D  bicubic  parametric  patch. 

(4)  Alignment  of  3D  surface  descriptions  based  on  sur¬ 
face  shape  using  heuristic  search  methods.  Genera¬ 
tion  of  numerical  models  of  objects  from  aligned 
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and  merged  surface  descriptions. 

(5)  Design  and  implementation  of  a  software  system  to 
accomplish  the  above  items  and  demonstration  of  its 
performance  using  several  modeling  examples. 

(6)  Design  and  implementation  of  a  software  system  for 
generation  of  synthetic  images  of  the  modeled  ob¬ 
jects?  expansion  of  the  pin-hole  camera  model 
traditionally  used  for  generation  of  synthetic  im¬ 
ages  to  include  several  parameters  of  an  actual 
optical  system  such  as  focus,  depth  of  field,  and 
motion  blur. 

7.3  Suggestions 

The  following  is  a  list  of  suggestions  for  possible 
continuation  and  expansion  of  the  work  presented  here: 

(1)  Segmentation  of  3D  numerical  object  models  into  3D 
surface  shape  and  volumetric  primitives  and  top- 
down  structuring  of  these  primitives. 

(2)  Matching  of  surface  segments  which  do  not  overlap, 
only  share  boundary  curves,  i.e.,  fragments  of 
broken  pottery. 

(3)  Detection  of  3D  surface  intersections  and  3D  sur¬ 
face  vertices  (points  where  three  or  more  surfaces 
intersect)  using  3D  curvature  search. 


(4)  Development  of  relational  matching  technique  of  3D 


objects  based  on  the  data  supplied  by  the  modeling 
process  described  here.  Automatic  conversion  of 
numerical  representation  into  symbolic  relations 
and  symbolic  surface  description. 

(5)  Expansion  of  the  optical  model  of  a  camera's  lens 
to  include  special-effect  filters  (such  as  star  anr’ 
diffraction}^  and  noise  caused  by  an  optical  sys 
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APPENDIX  A 


DEFINITION  OF  A  BICUBIC  PARAMETRIC  PATCH 
AND  A  SHEET  OF  PATCHES 

This  appendix  discusses  two  methods  of  constructing  bi¬ 
cubic  patches  from  the  surface  information  provided  by  the 
photogrammetric  technique  of  Chapter  2,  and  the  construction 
of  sheets  of  patches  from  individual  but  contiguous  patches. 
The  first  method  interpolates  a  patch  into  the  positional 
coordinates  of  four  control  points,  and  achieves  first- 
derivative  continuity  (C*)  across  patch  boundaries.  The 
second  method  constructs  a  patch  from  its  four  boundary 
curves;  its  continuity  across  patch  boundaries  depends  on 
the  continuity  of  the  boundary  curves  which  is  also  C1. 

A  bicubic  patch  is  a  cubic  polynomial  of  two  variables 
(parameters)  expressed,  in  matrix  form,  as: 

(A.  1 ) 

g(u,v)  »  £  u3  u2  u  l]]  bn  b12  bl3  bl4  v3 

b21  b22  b23  b24  y2 

b31  b32  b33  b34  v 

_  b41  b42  b43  b44  J  L 1  _ 

or,  for  short, 

(A. 2) 

g  (u,  v)  ■  U  B  Vb 

where  the  domain  of  the  parameters  ji,  and  y  is  the  unit 
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square:  0  £  u,  v  i  1.  In  equations  (A. 1-2)  g(u,v)  actually 
represents  a  component  of  the  vector: 


g  (u,  v) 


x  (u,  v) 
y  (u,  v) 
z  <u,  v) 


(A.  3) 


There  are  several  methods  of  computing  the  16  coeffi¬ 
cients  of  matrix  B  in  equation  (A.l),  such  as  the  Coons, 
Ferguson,  Bezier,  Hermite,  and  B-spline  methods 
[7,20,26,271.  The  Ferguson  method  is  useful  for  fitting  a 
patch  or  a  set  of  patches  into  existing  surface  data  points; 
the  Coons  method  fits  patches  into  boundary  curves,  whereas 
the  other  methods  are  more  useful  for  interactive  design  of 

i 

a  new  surface.  In  all  these  methods,  however,  B  has  to  be 
always  computed  from  16  items  of  information  about  the  sur¬ 
face  shape. 

A  Ferguson  patch  [271  is  constructed  from  surface  data 
of  four  control  points  as  follows;  let 

(A.  4) 

B  =  M  I  Mfc 


where 

”2-2  1  1 
-3  3  -2  -1 

M  * 

0  0  10 
10  0  0 


(A. 5) 
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is  a  matrix  of  coefficients  of  blending  functions,  and 

(A.  6) 

h  (0, 0)  h  (0, 1)  hv(0,0)  hv(0,l) 

h  (1, 0)  h  (1,1)  hv(l,0>  hv(l,l) 

hu(0,0)  hu(0,l)  huv(0,0>  huv(0,l) 
hu(l,0)  hu (1, 1)  huv(l,0)  huv (1,1) 

provides  the  description  of  the  surface  data  [Figure  A. 1(a)] 
with  the  following  short-hand  notation: 

dh  (u,  v) 

h  (u,v)  *  -  , 

du 

dh  (u,  v) 

h  (u,v)  =  -  ,  and 

dv 

d2h (u, v) 

h  (u,v)  »  -  , 

du  dv 

where  h(u,v)  is  the  positional  coordinate  of  a  control 
point,  hu(u,v)  is  the  surface  tangent  in  the  direction  of 
the  parametric  curves  y  -  constant,  hv(u,v)  is  the  surface 
tangent  in  the  direction  of  the  parametric  curves  .u  *= 
constant,  and  huv(u,v)  is  the  surface  cross-derivative  or 
twist.  The  values  of  hu(u,v),  hv(u,v),  and  huv(u,v)  can  be 
estimated  by  numerical-analysis  methods  [19]  from  the 
positional  data  h(u,v)  and  adopted  to  our  requirements  as 
follows: 
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h(u+ Au,v)  -  h(u-Aurv) 
2  Au 


,v 

.V 


hu<u,v)  = 


hv  (u,  v) 


h(u,v+Av)  -  h(u,v-Av) 
2  Av 


huv<u,v> 


[h  (u+  Au,  v+  Av) -h  (u-  Au,  v+  Av)  ]  - 
[h(a+Au,v-Av)-h(u-Au,v-Av)] 


4  Au  Av 


(A. 8) 


(A. 9) 


In  equations  (A. 7-8)  above/  the  surface  tangents  were 
computed  as  central  differences?  a  better  estimate  of  the 
tangents,  yielding  a  smoother  surface,  is  obtained  from  a 
weighted  average  in  a  3  x  3  point  neighborhood: 

(A. 10) 


hu(u,v) 


[h(u+Au,v+Av)+2h(u+Au,v)+h(u+Au,v-Av)J  - 
[h  (u-  Au,  v+  Av)+2h  (u-  Au,  v)+h  (u-  Au,  v-  Av)  ] 


8  Au 


hv(u,v) 


(A. 11) 

lh(u+ Au,v+ Av)+2h(u,v+ Av)+h(u-Au,v+ Av)l  - 
Ih(u+Au,v-Av)  +  2h(u,v-Av)  +  h(u-Au,v-  Av)  ] 


8  Av 


In  summary,  a  Ferguson  patch  is  completely  defined  by  [ 
h(u,v)  hu(u,v)  hv(u,v)  huv(u,v)  ]  evaluated  at  each  of  the 
four  corner  control  points  of  the  patch  and  by  the  blending 
coefficients  of  (A. 5).  For  two  adjacent  patches  to  have 
positional  and  derivative  continuities,  they  must  share  the 


surface  data  defined  at  their  common  control  points. 

A  Coons  patch  [20]  is  constructed  from  four  boundary 
curves  c(u,0),  c(u,l),  c(0,v),  and  c(l,v)  [Figure  A. 1(a)]  by 
linear  interpolation  between  the  opposite  pairs  of 
boundaries : 

g  (uf  v)  =  £  (1-u) 

+  £  (1-v) 

-  £  (1-u) 

The  four  boundary  curves  in  equation  (A. 12)  can  be 
computed  by  the  reconstruction  method  given  in  Appendix  B. 
Those  curves  have  C*  continuities  in  their  2D  projections  as 
well  as  in  the  3D  reconstructions.  Therefore,  a  patch 
constructed  from  such  boundary  curves  will  also  have  C1 
continuity  across  its  boundaries. 

A  sheet  of  composite  bicubic  patches  consists  of  a  set 
of  contiguous  patches  with  at  least  C*  continuity  across 
their  boundaries.  The  patches  in  a  sheet  are  parametrized 
by  a  single  parameter  coordinate  system  P(u,v,w).  This  co¬ 
ordinate  system  [Figure  A. 1(b)]  contains  two  sets  of  orthog¬ 
onal  lines  which  define  an  integer  square  grid.  These 


»D 


1 


ul 


(A. 12) 


c  (0,  v) 
c  (1,  v) 

c  (u,  0) 
c  (u,  1) 

C (Or  0)  C (0, 1) 
c  (1,0)  c(l,l) 


(1-v) 

V 


tat 
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lines,  mapped  on  3D  surfaces,  become  the  patch  boundary 
curves  and  their  intersections  become  patch  control  points. 
The  lines  define  an  absolute  parameter  coordinate  system  for 
0  £  ji  i  fi,  and  0  i  .2  i  JS,  while  within  a  single  patch  there 
is  a  relative  parameter  coordinate  system  for  0  <  u,  v  <.  1. 
Using  the  notation  of  the  absolute  parameter  coordinate  sys¬ 
tem,  a  patch  is  defined  by  the  four  control  points 
h(m,n),  h(m+l,n),  h(m,n+l),  and  h  (m+1 ,  n+1) ,  or  by  the  four 
boundary  curves  c(m,v),  c(m+l,v),  c(u,n),  and  c(u,n+l)  for  m 
<l  11  1  m+1,  and  n  i  z  <  u+1. 


APPENDIX  B 


RECONSTRUCTION  OF  A  PARAMETRIC  CURVE 
FROM  MULTIPLE  PROJECTIONS 

In  Section  2.3  we  reconstructed  3D  surface  points  from 
multiple  2D  perspective  projections.  These  points  were  then 
used  as  the  control  points  for  parametric  patches.  In  this 
appendix  we  develop  an  analogous  method  which  reconstructs  a 
3D  parametric  space  curve  from  its  multiple  2D  projections 
[Figure  B.U.  Such  a  curve  can  directly  become  one  of  the 
four  boundary  curves  of  a  surface  patch. 

A  3D  cubic  curve  in  an  object  coordinate  system 
0(x,y,z,l)k  is  specified  by: 

c(w)k  = 

or,  for  short, 

(B.2) 

c(w)R  «  Ck  W 

Similarly,  a  2D  projection  of  this  curve  into  an  image 
coordinate  system  Q(r,s,l)^  is: 


x  (w) 

C11 

c12 

c13 

1 - 

H 

O 

y  (w) 

s 

C21 

c22 

c23 

c24 

z  (w) 

k 

_  C31 

c32 

c33 

c34  _ 

[•] 


Lsiw;j.  Lc21  c22  c23  c24ji  w 


or,  foe  short. 


clv)^  =  W 


(B.4) 


The  mapping  from  the  object  coordinate  system  0k  to  the 
image  coordinate  system  Qi  is  T{Ok-»Qi).  Each  perspective 
projection  i  of  the  2D  curve  has  two  components:  r(w)^,  and 
s(w)jl;  the  3D  curve  has  three  components:  xtw)^,  y(w)k, 
and  A  component  of  a  cubic  curve  is  constructed  from 
four  items  of  curve  information.  These  are,  typically,  the 
position  coordinates  of  the  two  end  points  (h(0),h(l))  and 
the  slope  of  the  curve  at  these  points  (t^  (0) ,1^  (1) ) .  Here, 
h(w)  refers  to  the  positional  coordinate  of  the  curve  at  w, 
and  hw  (w)  refers  to  the  slope  of  the  curve  at  The  param¬ 
eter  w  is  defined  for  0  i  a  i  1. 


A  Ferguson  curve,  similar  to  a  Ferguson  patch  of 
Appendix  A,  is  constructed  from  the  curve  information,  and  a 
matrix  of  coefficients  of  blending  functions: 


ii 


c  (w) 


[  h (0)  h (1)  hw(0)  hw(l>  ] 


(B.5) 


"2-3  0  1 ~ 

w3 

-2  3  0  0 

w2 

1-210 

w 

1  -i  0  0  _ 

_  1 

It  is  necessary  to  point  out  that  the  definition  of  the 
parameter  a  must  be  consistent  in  all  projections  i;  that 
is,  a  value  of  a  yields  2D  projections  of  the  same  point  in 
3D.  The  system  of  equations  (2.3)  can  be  now  rewritten  by 
the  substitution  of  x(w)k,  y(w)k,  z(w)k,  r(w)if  and  sfw^ 


Xk,  yk. 

zk,  r,. 

and  s^. 

respectively  ( 

.  as: 

(w) 

di2(w) 

di3(w) 

d14 (w) 

~x(w) 

2i  (w) 

d22(w) 

d24(w) 

d24 _ 

y  (w) 

z  (w) 

(B.  6) 


=  0 


where  d^p(w)  and  d2p(w)  are  defined  by: 


dip(w) 


fc3p  (  clliw  +  c12iw  +  c13iw  +  c14i 


)  - 


(B.  7) 
fclp 


d2p(w)  18  fc3p  (  c21iw3  +  c22iw2  +  c23iw  +  c24i  } 


-  t 


2p 


for  p  ®  1,2, 3, 4. 

Also,  all  £'s  and  I's  in  equation  (B.7)  apply  to  trans- 
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and  z (w)  are  defined  in 


formation  i,  and  x(w),  y(w), 

equation  (B.l).  We  need  to  solve  the  system  of  equations  in 
(B. 6)  for  the  twelve  coefficients  of  given  in  (B.l).  The 
matrix  elements  of  equation  (B.7)  and  the  vector  elements  of 
equation  (B.l)  are  multiplied  in  equation  (B.6)  and  the  co¬ 
efficients  of  the  curve  (B.l)  are  factored  out  to  give  the 
following  system  of  two  linear  equations  with  twelve 
unknowns: 

(B.8) 


d11(w)w3  dn(w)w2 
d2i (w) w^  d22 (w)  w2 


d^3 (w) w 
d23 (w) w 


d13 

d23 


d14 

d24 


(w) 

(w) 


Each  of  the  two  equations  in  (B.8)  can  be  expanded  into 
four  linearly  independent  equations  by  substituting  four 
linearly  independent  vectors  W  (values  wlf  &2,  #3*  and  j*4) 
for  the  parameter  w,  because  the  cubic  basis  vector  w  yields 
four  linearly  independent  vectors.  Therefore,  for  a  single 
projection  i.  there  is  a  system  of  eight  equations  with 
twelve  unknowns,  expressed  as: 


06 


d13(wl} 


d14 (w3) 


dll(wl)wl 


dll(w2)w2 


di3(w2) 


21(w3)wi  •** 

d23(w3} 

21 ^w4  ^ w4  • • • 

d23(w4) 

C11 

c12 


c33 

L  c34  J 


d14 (w2) 


d 24  (w3> 


d24(w4}  J 


or,  for  short/ 


(B. 10) 


DC'  =  D' 

In  general,  for  12  2  projections/  there  j.s  a  system  of 
8 • I  equations  of  twelve  unknowns: 

(B.ll) 


1 

ol 

M 

_ 1 

1 

Ol 

_ 1 

Di 

C  * 

D'i 

1 - 

Ol 

M 

1 

Ol 

M 

1 

System  (B.ll)  is  solved  by  the  least-squares  method 
which  yields: 

(B.12) 


- 1 

Ol 

M 

_ 1 

t 

l 

ol 

M 

_ 1 

t 

1 

ol 

_ 1 

Di 

C'  = 

Di 

D'i 

1 - 

Ol 

M 

_Di_ 

_Dl_ 

_D’l_ 

This  curve  reconstruction  method  is  particularly  useful 
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when  the  2D  projections  of  the  curve  are  computed  from  in¬ 
formation  other  then  four  position  points,  i.e.,  two  end 
points  and  two  slopes  as  in  equation  (B.5),  or  when  the 
points  in  one  projection  cannot  be  matched  with  their  corre¬ 
sponding  points  in  the  other  projections.  This  method  can 
be  generalized  to  parametric  curves  with  an  arbitrary-order 
basis  vector  W. 
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APPENDIX  C 


SURFACE  REPRESENTATIONS  AND  OPERATIONS 


This  appendix  contains  a  summary  of  representations  of 
the  four  types  of  3D  surface  elements  (planar/  spherical/ 
quadric/  and  bicubic)  and  the  five  basic  geometric 
operations  (3D  rigid  transformation,  surface-normal  vector, 
surface-normal  curvature,  intersection  with  3D  line,  and 
surface  existence)  on  them.  This  is  followed  by  individual 
descriptions  of  the  representations  and  the  operations  for 
each  type  of  surface  elements. 

C.l  Summary 


C. 1.1  Representations  of  Surface  Elements 


(a)  algebraic: 


(C.l) 

f(x,y,z,w)  *  0 

bounds:  other  algebraic  representations 

arranged  in  a  boolean  tree 


(C.  2 ) 

(b)  parametric:  g(u,v)  *  [x(u,v),  y(u,v),  z(u,v),  w(u,v)^| 

bounds:  limits  on  parameters  ji  and  y 
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C.1.2  Geometric  Operations 


(a)  transformation  from  object  coordinate  space  0(x,y,z,w)k 
to  object  coordinate  space  0(x,y,z,w)k,  by  T{Ok-»Oki}: 

(C.3) 


T  (0k— »  }  = 


fcll 

fc12 

fc13 

rr 

»-■ 

fc21 

fc22 

fc23 

fc24 

fc31 

fc32 

fc33 

fc34 

0 

0 

0 

1 

The  transformation  T(Ok-»Oki}  is  a  function  of  the 
following  transformation  parameters: 

(C.4) 


translations: 

tx. 

ty,  tz 

rotations: 

rx, 

ry,  rz 

translations  of 

rotations: 

trx, 

try. 

trz 

order  of  rotati 

ons: 

irx, 

iry. 

irz 

scale  factors: 

sx. 

sy,  sz 

translations  of 

scale  factors: 

tsx, 

tsy. 

tsz 

The  transformation  parameters  are  normally  written  in 
matrix  form.  The  translation  matrix  is: 


f  <tx,  ty,  tz)  = 


1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

1 

0 


-tx 

-ty 

-tz 

1 


(C.  5) 


the  rotation  matrices,  in  the  irz.  i rv.  and  i rx  order  of 


2 


[i 


precedence  are: 


R (rx, ry, rz) 


I 


1  0 

0 

0  cos(rx) 

sin  (rx) 

0  -sin(rx) 

cos  (rx) 

0  0 

0 

cos  (ry) 

0 

-sin  (ry) 

0 

1 

0 

sin  (ry) 

0 

cos (ry) 

0 

0 

0 

cos  (rz) 

sin(rz)  0 

-sin  (rz) 

cos(rz)  0 

0 

0  1 

0 

0  0 

and  the  scale  matrix  is: 

r  sx 


i 


S (sx, sy, sz) 


0 

0 

0 


0 

sy 

0 

0 


0 

0 

sz 

0 


0 

0 

0 

1 


(C.6) 


(C.7) 


The  complete  transformation  is  then  a  concatination  of 
scale  transformations  around  point  R (tsx, tsy , tsz) , 
followed  by  z,y,x  rotations  around  point  R (trx, try, trz) , 
followed  by  translation  to  point  h(tx,ty,tz): 

(C.8 ) 


T(Ok-»Oki}  -  T(tx,ty,tz) 

T (-trx,-try,- trz) R (rx, ry, rz) T (trx, try, trz) 
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The  inverse  matrix  T(Ok— »ok « transforms  from  object 
coordinate  space  0<x,y,z,w)k.  to  object  coordinate  space 
0 (x, y, z, w) k,  that  is: 

(C.9) 

T{0k->0ki}“1  =  T{0ki-»0k> 

(b)  surface-normal  vector  at  point  "K(x,y,z)  of  algebraic 
representation  f(x,y,z): 

(C.10) 

N  *  £df/dx  df/dy  df/dz^] 

surface-normal  vector  at  point  Ti(u,v)  of  parametric  rep¬ 
resentation  g(ufv): 

(C.ll) 

N  =  g  (u,  v)  u  X  g  (uf  v)  y 

■  £ dx/du  dy/du  dz/du]]  X  [[  dx/dv  dy/dv  dz/dv^| 


(c)  normal  surface  curvature  g,  at  point  Ti(u,v)  of  parametric 
representation  g(u,v)  in  the  direction  of  the  surface 
curve  g(u(t),v(t))  is  given  by: 

(C.12) 

D^U^  +  2D2UV  +  DjV^ 

q  3  - 

E^U^  +  2E2UV  +  EgV^ 


where 


(C. 16) 


algebraic: 


parametric: 


transformation: 


f (x, y, 2)  =  a0x  +  axy  +  a2z  +  a3  =  0 

-  C«0  »2  Ml  ["*' 

Y 
z 
1 

bounds:  lists  of  edges  and  vertices  [Figure 
C. 1 (a) ] 

(C. 17) 


g  <u, v)  =  C  b0  bl  ^2  D 


u 
v 
1 

bounds:  limits  on  parameters  u  and  v, 

possibly  combined  with  lists  of  edges  and 
vertices 

(C. 18 ) 


1 

01 

0 

1 - 

0 

m 

al' 

=  T (0— ^0 ' ) 

al 

a2  ' 

a2 

_  a3 '  _ 

_a3_ 

(the  algebraic  representation  in  equation 
(C.16)  is  pre-multiplied  by  T{0-*0'}) 


Figure  C.l 


Surface  representation:  (a)  planar  face 


(C.  1 9) 


normal  vector:  ^  al  a2  H 


N  =  Cb0x  b0y  b0zZI  *  [blx  bly  blzD 


(C. 20) 


intersection:  f (lv (t) ,1„ (t) ,1_ (t) )  =  0 


(obtained  by  substitution  of  parametric  line 
representation  (C.13)  into  algebraic  surface 
representation  (C.16),  a  linear  equation 
with  up  to  one  possible  solution  [Figure 
C.l(a)];  the  number  of  valid  intersections 
of  a  line  from  the  point  of  intersection  to 
infinity  in  the  plane  with  all  edges  of  the 
face  then  determines  whether  the  inter¬ 
section  point  is  inside  oc  outside  the 
planar  face:  if  the  numbei  is  odd  the  point 
is  insidef  otherwise  it  is  outside;  this 
method  allows  a  planar  fact  bounded  by  a 
closed  list  of  edges  to  contain  any  number 
of  not-intersecting  holes  or  protrusions, 
each  similarly  bounded  by  a  closed  list  of 


(C.21) 


algebraic: 


parametric: 


transformation: 


normal  vector: 


f(x,y,z)  =  (x-xc)2+ (y-yc)2+ (z-zc)2~r2  *  0 
bounds:  none  -  the  representation  is  self- 
bounding 


(C.22) 

x(u»v)  =  xc  +  r  cos  (u )  cos  ( v) 
y(u,v)  =  yc  +  r  cos(u)  sin(v) 

2 (u, v)  *  zc  +  r  sin  (u) 
bounds:  -pi/2  £  ^  £  pi/2 
0  i  2  si 

(C. 23) 


"v~ 

~xc~ 

yc' 

»  T{0->0' } 

yc 

zc' 

2C 

_  wc '  _ 

_  wc_ 

(the  center  point  Ti  (xc,  yc,  zc,  wc)  of  the 
sphere  is  pre-multi plied  by  T(0-»0'}) 


N  =  [x-xc  y-yc  z-zcn 
or 


(C. 24) 


2 


f 

intersection: 


-sin(u)  cos  (v) 

-cos(u)  sin(v) 

N  ■ 

-sin(u)  sin(v) 

X 

cos(u)  cos  (v) 

cos  (u) 

0 

(C.  25) 

f  (lx(t),ly(t)flz(t))  =  0 

(obtained  by  substitution  of  parametric  line 
representation  (C.13)  into  algebraic  surface 
representation  (C.21),  a  quadric  equation 
with  up  to  two  possible  solutions  (Figure 
C.l(b)l;  if  the  sphere  serves  as  a  bounding 
volume  it  is  not  necessary  to  compute  the 
exact  intersection  points,  therefore  only 
the  line  to  center-of-sphere  distance  is 
computed  and  compared  with  the  radius  of  the 
sphere) 


C.4  Quadric  Surface 


(C.26) 


algebraic: 


f(x,y,z)  = 


2  2  2 

a8x*  +  a^y*  +  a2z*  + 

a3xy  +  a4xz  +  a5yz  + 

agx  +  a-jY  +  ag  z  +  ag  = 


=  |x  y  z  1J  a8  83/2  a4/2  ag  x 

a^/2  a^  ag/2  a-j  y 

a4/2  ag/2  a2  a8  z 

a6  a  7  a8  a9  1 

“  [x  y  Z  1]  A  [x  y  z  ljt 
bounds:  boolean  tree  of  quadric  surfaces 


parametric: 


(C.27) 


g(u,v)  «  [u2  u  l]  To 


b22  b23 


b31  b32  b33 


*  U  B  ^ 


bounds:  limits  on  parameters  u  and  y 


(C. 28 ) 


transformation:  A' 


T(0— X?'  )_1  A  (T (0— >>0 • } “1 )  fc 


(matrix  A  of  equation  (C.26)  is  pre¬ 
multiplied  by  the  inverse  of  T{0-»0'>  and 
post-multiplied  by  the  transpose  of  the 


normal  vector: 


intersection: 


inverse  of  T(0-»0'}  [41]) 


N 


2agx  +  a3y  +  a4z  +  ag 
2a3y  +  a3x  +  a5z  +  a7 
2a2z  +  a4x  +  a5y  +  ag 


or 


(C.  2  9) 


N  = 

r  i 

(2bi3U+b22v+b23)x 
(2b33U+b22v+b23)  y 

X 

(2b31v+b22u+b32)x 

(2b31v+b22 u+b32  ^  y 

(2b13u+b22v+b23)  z 

(2b31v+b22u+b32)  z 

(C. 30) 


f  <lx(t),ly(t),lz(t))  =  0 

(obtained  by  substitution  of  parametric  line 
representation  (C.13)  into  algebraic  surface 
representation  (C.26),  a  quadric  equation 
with  up  to  two  possible  solutions  [Figure 
C.2  (a)] ) 


c  Pa 


algebraic: 


none 


parametric: 


g  <u,  v) 


*  £ u^  u  1 J 


bll  b12  b13  b14 
b21  b22  b23  b24 
b31  b32  b33  b34 
b41  b42  b43  b44 


(C. 31) 

ir..3i 


■  0  I  vb 

bounds:  0  1  u,  v  i  1  in  relative  parameter 

coordinates 

(C.  32 ) 

transformation:  gx(u,v)'  =»  T(0— »0'>  gx(u,v) 

gy(u,v)'  =  T{0-»0'}  gy(u,v) 

gz(u,v)'  =  T{0-»0'}  gz(u,v) 

(matrix  B  of  each  component  in  equation 
(C.31)  is  pre-mul tiplied  by  T(0-»0'}) 


normal  vector: 


intersection: 


(C. 33) 

N  =  g(u,v)u  X  g(u,v)v 

(C. 34) 

g  (u,  v)  -  T  (t)  *  0 

or  when  written  as  three  simultaneous  cubic 
equations: 


(C. 35) 


gx(u#v)  -  ix(t)  =  o 

gy(u,v)  -  ly(t)  ■  0 

gz(u,v)  -  lz(t)  =  0 

This  system  of  non-linear  equations  is 

numerically  solved  by  a  modified  Newton- 
Raphson  method  with  tests  for  equation  con¬ 
vergence  to  the  nearest  intersection  solu¬ 
tion#  and  constraints,  limiting  the  parame¬ 
ters  to  within  the  patch  and  the  bounding 
volume  [Figure  C.2(b)].  Intersection  of  the 
line  with  a  planar  approximation  of  the 
patch  or  its  subpatch  provides  the  initial 
estimate  of  the  three  unknowns. 

The  solution  proceeds  as  follows;  let: 

(C. 36) 

E(t»u,v)  ■  gx(u,v)  -  lx(t) 

F(t,u,v)  *  gy  (u,  v)  -  ly(t) 

G(t,u,v)  -  gz(u,v)  -  lz(t) 


Then  by  first-order  Taylor  expansion  of  the 
system  of  equations  in  (C.36)  we  can  incre- 
ment'.iy  solve  for  the  parameters  t,  u,  and 


(C. 37) 

E  dE/du  dE/dv 
F  dF/du  dF/dv 
G  dG/du  dG/dv 


I  D  (tf 

u. 

v)  | 

( 

dE/dt 

E 

dE/dv 

dF/dt 

F 

dF/dv 

dG/dt 

G 

dG/dv 

|  D  (t 

,u,v)  | 

(C. 

.39) 

dE/dt 

dE/du 

E 

dF/dt 

dF/du 

F 

dG/dt 

dG/du 

G 

I  D(t,u,v)  | 
where  |D(t,UrV)|  is  the  determinant  of  the 
Jacobian  matrix  defined  by: 


D  (t,  u,  v) 


An  outline  of  the  algorithm  follows: 


(C. 41 ) 


)rocedure  csystem (t, u, v) ; 

keqin 

I  test  magnitudes  of  the  three  equations] 


while  IE(t,u,v)l  >  epsi  I 
I F (t, a, v) I  >  eps^  I 
I G (t, u, v) t  >  eps2  do 
begin 

[test  for  vanishing  determinant  of 

Jacobian  matrix] 

dd  :«  determinant (D  (t, u, v) ) ; 

if  Idd I  <  eps2  then  return (fail) ? 

dt  :=  [from  equation  (C. 37 ) 1 y 

du  :*  [from  equation  (C.38)l; 

dv  :=  [from  equation  (C. 39)]; 

[test  boundary  condition,  if  outside 
clip  to  boundary  values,  if  previously 
on  boundary  then  fail] 
if  (t-dt  >  t„av)  I  (t-dt  <  I 


(t-dt 
(u-du 
(v-dv 
begin 
if  ct 
(u 
(v 


kmax} 

“max 

vmax> 

*ininj 

“min 


(u-du  <  ul 
(v-dv  <  vm 

•  W 

“  -  “max 
(v  =  Vmav) 


'-min^  j 
Umin! 

vmin)  then 


(v  -  v-i")  I  (v  =  v-  )  iten 

return  (fail) : 

else  clip  to  boundary  values; 
[test  convergence  to  the  nearest 
solution,  if  increments  become  too 
small  exit  with  failure] 
while 

I E  (t-dt, u-du, v-dv) I > I E (t, u, v) I  I 
I F (t-dt, u-du, v-dv) I > I F  (t, u, v) I  I 
lG(t-dt, u-du, v-dv)  l>IG(t,u,v)  I 

is 

baain 

dt  :=  dt  /  2; 

du  :■  du  /  2; 

dv  :=  dv  /  2; 

[test  if  at  local  minimum  or 
maxim™] 

1  f  J..2  /  Ane  ^ 


local  minimum 


end; 

end? 

[compute  new 
tion] 

t  :*  t  -  dt? 
u  ;■  u  -  du? 
v  :*  v  -  dv? 
end? 

‘turn (t, u, v) ? 


if  dt^+du^+dv^  <  epsj 

then  return  (fail ) ? 


increment 


the  solu- 


APPENDIX  D 


SOFTWARE  COMMANDS 


This  appendix  contains  descriptions  of  the  available 
commands  in  the  four  software  command  processors,  SURE, 
SUED,  ALIGN,  and  STRAW,  described  in  Chapter  5.  The  syntax 
of  the  commands  is  given  by: 

(D.l) 

COMMAND:  <parameter_list>; 

where 

<parameter_list>  ::=  <parameter_list>  <operator> 

<parameter>  I  <parameter> 
<parameter>  <identifier>  I  <constant>  I 

(keyword) 

<operator >  '  '  I  I  '  '  I  'I'  I  '+'  I 

I 

Here,  <identifier>  refers  to  the  name  of  data  in  a  proc¬ 
essor's  data  base,  or  to  the  name  of  a  directory  or  a  file; 
<constant>  refers  to  a  numerical  constant  in  either  integer 
or  floating-point  format;  and  (keyword)  contains  a 
command's  keyword. 

D.l  Processor  SURE 


TAPEIN:  <directory_name>,  <8-bit_file_name>. 
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<16-bi  t_file_name>; 

reads  a  digitized  image  from  the  current  position  of 
a  magnetic  tape  and  creates  an  8  bit/pixel  and, 
optionally,  a  16  bit/pixel  image  file. 

HISTOGRAM:  <directory_name>,  <file_name>; 

computes  a  histogram  of  a  16  bit/pixel,  8  bit/pixel 
or  1  bit/pixel  (binary)  image. 

CONVERT:  <di rectory_name>,  <in_f ile_name>,  <out_file_name>? 

converts  either  a  16  bit/pixel  image  to  an  8 
bit/pixel  image,  or  an  8  bit/pixel  image  to  a  binary 
image. 

WINDOW:  <directory_name>,  <in_file_name>,  <out_file_name>; 
creates  a  window  image  from  an  input  image. 

FFT2D:  <di rector y_name>,  <inufile_name>,  <out_file_name>; 

computes  the  2D  fast  Fourier  transform  of  an  8  or  16 
bit/pixel  image. 

FILTER:  <f il ter_parameters> ; 

filters  an  image  in  the  frequency  domain. 

IFFT2D:  <directory_name>,  <in_file_name>,  <out_f ile_name>; 

computes  the  2D  inverse  fast  Fourier  transform  and 
generates  an  8  or  16  bit/pixel  image. 

MFFT2D:  <di rectory_name>,  <in_f ile_name>,  <out_file_name>; 

converts  a  2D  Fourier  transform  of  an  image  into  a  8 
bit/pixel  image  of  logarithmically  scaled  magnitudes 
of  the  Fourier  transform. 

CLEAN:  <di rectory_name>,  <in_file_name>,  <out_file_name>; 
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smooths  either  an  8  bit/pixel  or  a  binary  image  with 
a3  x  3  or  a5x5  pixel  operator  in  the  spatial 
domain. 

PLOT:  <directory_name>,  <in_file_name>; 

plots  an  8  bit/pixel  or  a  binary  image  on  the 
electrostatic  printer-plotter.  The  8  bit/pixel  im¬ 
age  is  dithered. 

COPY:  <directory_name>,  <in_file_name>; 

copies  a  window  of  an  8  bit/pixel  or  a  binary  image 
to  the  frame  buffer  of  the  raster  terminal. 

LOAD:  <identif ier>,  <directory_name>,  <file_name>; 

loads  an  image  or  partially  processed  data  into 
SURE. 

EXTRACT:  <identifier>; 

extracts  the  image  projections  of  the  corners  of 
camera  calibration  marks  from  a  binary  image. 
MATCH_MARKS:  Cidentif ier >; 

matches  the  2D  image  projections  of  the  corners  of 
camera  calibration  marks  with  their  3D  coordinates. 
TRANSFORMATION:  <identif ier>; 

computes  the  image  transformation  matrix  by  matching 
a  3D  model  of  the  calibration  marks  with  their  2D 
projections  in  the  specified  image. 

THIN:  <identif  ier>; 

iteratively  thins  a  binary  image  until  a  line  pixel 
has  at  most  only  two  neighbors. 


RECONSTRUCT:  <image_l>,  <image_2>,  ....  ; 

reconstructs  a  3D  surface  representation  from  infor¬ 
mation  obtained  from  two  or  more  images. 

SAVE:  <identif ier>,  <directory_name>,  <file_name>; 

saves  an  image  or  partially  processed  data  from  SURE 
to  a  file. 

LIST:  (marks  I  grid  I  nodes  I  curves} ; 

lists  contents  of  the  specified  data. 

DISPLAY:  <identif ier>; 

displays  the  specified  data  on  the  vector-drawing 
graphics  terminal. 

D.2  Processor  SUED 

SPHERE:  <identif ier >,  <color>,  <x>,  <y>,  <z> ,  <radius>; 

defines  the  name,  color,  location  and  size  of  a 
sphere. 

VERTEX:  Cidentif ier >,  <x>,  <y>,  <z>; 

defines  a  3D  point  that  can  be  used  either  as  a 
vertex  point  in  planar  faces  or  as  a  control  point 
in  bicubic  patches. 

PLANE:  <identif ier >,  <color>,  <vertex_list>  &  <vertex_list> 

&  ....  i 

defines  a  planar  face  that  is  bounded  by  a  list  of 
vertices  specifing  outer  edges,  and,  optionally,  one 
or  more  lists  of  vertices  specifing  inner  edges 


(holes).  A  bounding  volume  is  cast  around  the 
planar  face. 

POLYHEDRON:  <identif ier>,  Cplane>,  <plane>  ....  ; 

defines  a  set  of  planar  faces  which  form  a 
polyhedron-like  object  (holes  in  faces  are  allowed 
and  the  object  does  not  have  to  be  closed).  A 
bounding  volume  is  cast  around  the  object. 

PATCH:  Cidentif ier >,  Ccolor>,  <list_of_control_points>, 
<r_min>; 

defines  a  bicubic  patch  that  is  fitted  into  16  3D 
control  points  which  are  defined  by  the  VERTEX 
command.  Such  a  patch  will  usually  have  at  most  C° 
continuity  with  any  adjacent  patches. 

QUADRIC_COEFFS :  Cidentif ier >,  Ccolor>,  (actual  I  auxiliary) , 

<aQ>,  <a^>,  ....  Ca^^r 

defines  the  name,  color,  type  and  algebraic  repre¬ 
sentation  of  a  quadric  surface.  The  type  of  the 
surface  may  be  auxiliary  -  used  only  as  a  bounding 
surface  -  or  normal. 

QUADRIC_BOUNDS :  Cidentif ier>,  Cbounding_tree>; 

defines  a  single  level  AND-OR  logical  tree  of 
bounding  quadric  surfaces. 

COLOR:  cidentif ier>,  Cr_ambient>r  Cg_ambient>,  cb_ambient>, 
Cambient_coef f >,  Cdif f use_coef f >, 

Creflection^coef f >,  ctransmission_coef f >, 

Cref raction_coef f >,  Cspecular_coef f >, 


cglossiness_coef f >; 

defines  the  ambient  red,  green  and  blue  intensities, 
and  the  ambient,  diffuse,  reflection,  transmission, 
refraction,  specular  and  glossiness  coefficients  of 
a  color. 

LOAD:  Cidentif ier>,  cdi rectory_name>,  <file_name>; 

loads  a  surface  representation  from  a  file  into  the 
data  base  of  the  processor. 

SAVE:  Cidentif ier >,  <di rectory_name>,  <file_name>; 

stores  a  surface  representation  into  a  file  from  the 
data  base  of  the  processor. 

SYMMETRY:  Cidentif ier >,  Cfrom_x>,  Cfrom_y>,  Cfrom_z>, 
cto_x>,  cto__y>,  cto_z>,  Cp_increment>, 
Cc_increment>; 

attempts  to  convert  a  sheet  of  bicubic  patches  into 
a  generalized  cylinder-like  shape. 

QUADTREE:  cidentif ier_guadtree>,  Cidentifier_patches>; 

converts  a  sequential  list  of  bicubic  patches  into  a 
structured  quadtree  of  the  patches  with  bounding 
volumes  and  normal  vectors. 

SET-DISPLAY:  (initialize),  (blink  I  noblink > ,  (fast  I  normal ) , 
(solid (dotted  I  dashed) ,  ( came ra_st and  I maxim urn) ; 
defines  the  display  mode  of  the  subsequently  dis¬ 
played  data. 

DISPLAY:  (earner a_stand I  nodes  I  patches (quadtrees  I  cylinders) ; 

displays  the  specified  part  of  the  currently  defined 
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data 


STATUS_DISPLAY; 

lists  the  current  status  of  the  display  device. 

LIST:  {nodes  I  patches [quadtrees Icylindets} ; 

lists  the  specified  part  of  the  currently  defined 
data. 

VERIFY :  (nodes  I  patches Iquadtrees I  cylinders} ; 

verifies  the  integrity  of  the  specified  data. 

MODIFY:  <node>,  <x>,  <y>,  <z>,  <pointer_list>; 

adds  a  new  node  or  modifies  the  contents  of  an 
existing  node  of  a  parametric  sheet. 

SMOOTH:  <sheet>,  <u_order>,  <v_order>; 

smooths  a  sheet  of  bicubic  patches  with  an  orthogo¬ 
nal  bivariate  polynomial  of  the  specified  orders. 

TRANSFORM:  <identif ier>,  <tx>,  <ty>,  <tz>,  <rx>,  <ry>,  ....; 

creates  a  3D  transformation  matrix  composed  from 
translations,  rotations  in  arbitrary  order  around  an 
arbitrary  point,  and  scale  changes  around  an 
arbitrary  point. 

RELATION:  <identif ier >.,  <parameters>; 

creates  a  relational-node  entry. 

BOUND:  <identifier>,  <parameters>; 

creates  a  boundary-volume  entry. 

INSERT:  Cidentif ier >,  <node>,  <transform>; 

inserts  a  graph  node  into  the  graph  data  structure. 

INITIALIZE;  STATUS;  HELP;  EXIT; 
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are  the  obvious  and  self-explanatory  utility 
commands. 

D.3  Processor  ALIGN 

LOAD_SURFACE :  <identif ier >,  <di rectory_name>,  <file_name>; 

loads  the  first  or  the  second  structured  surface 
into  the  data  base. 

SET_TRANSFORM:  {txlno_tx},  <max_tx>,  <min_tx>,  ....  ; 

defines  transformation  parameters  which  can  be 
modified  during  the  search  process  and  sets  limits 
of  their  ranges. 

MODE:  (debug  I  nobug} ,  (display  I  noplay} ; 

enables  and  disables  debug  and  display  modes  of  the 
search  process. 

MATCH ; 

aligns  the  two  currently  defined  surfaces  into  a 
common  object  coordinate  system. 

MERGE; 

merges  the  two  aligned  surfaces,  defined  in  the  data 
base,  into  a  common  parameter  coordinate  space. 

SAVELSURFACE :  <di rectory_name>,  <file_name>; 

saves  the  aligned  and  merged  surfaces  into  a  file. 

LOAD_SEARCH:  <di rectory_name> ,  <file_name>; 

loads  the  search  data  base  generated  by  the  ALIGN 
command  from  a  file. 
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SAVE_SEARCH :  <di rector y_name>,  <file_name>; 

saves  the  search  data  base  to  a  file. 

TRACE:  <search_node>; 

traverses  the  search  tree  of  the  alignment  process 
to  the  specified  node  and  displays  the  surface 
transformation  of  each  node  in  the  path. 

TRACE_GOAL ? 

traverses  the  search  tree  of  the  alignment  process 
to  the  goal  node  and  displays  the  surface  transfor¬ 
mation  of  each  node  in  the  path. 

REPLAY? 

repeats  a  previously  executed  alignment  process  by 
traversing  the  nodes  of  its  search  tree  in  the  order 
in  which  they  were  generated. 

INITALIZE;  STATUS;  HELP?  EXIT; 

are  again  the  obvious,  indispensable,  and  self- 
explanatory  utility  commands. 

D. 4  Processor  STRAW 

CAMERA:  <x>,  <y>,  <z>,  <rot_x>,  <rot_y>,  <rot_z>, 

<focal_length>,  (perspective lorthographic} ; 
defines  the  position  of  the  camera  model  (center  of 
projection)  in  the  object  coordinate  system,  rota¬ 
tions  of  the  optical  axis  of  the  camera  model,  the 
focal  length  of  its  lens,  and  the  type  of  pro- 


jection.  In  the  orthographic  projection  the  focal 
length  only  scales  the  rays'  parameter;  all  rays 
are  parallel  to  each  other  and  perpendicular  to  the 
image  plane. 

FRAME:  <r_size>,  <s_size>,  <r_pixels>,  <s_lines>, 

<i_pixel>f  <f_pixel>,  <i_line>,  <f_line>; 
defines  the  size  of  the  image  frame  in  object  coor¬ 
dinate  system  units,  in  pixel  units,  and  defines 
window  in  the  image  to  be  actually  computed.  The 
size  of  the  image  plane  and  the  focal  length  of  the 
camera  model  together  define  the  viewing  pyramid. 

SAMPLE:  <levels>,  <subdivision_coef f s>; 

specifies  a  maximum  number  of  levels  that  a  pixel 
can  be  subdivided  to  perform  image  antialiasing.  A 
subdivision  coefficient  for  each  level  specifies  the 
amount  by  which  the  average  intensity  of  a  pixel 
area  can  differ  from  any  one  of  its  samples  before 
the  area  is  subdivided. 

SHADE:  <max_ray>,  <min_hide>,  <mirushadow>,  <max_ref lect>, 
<roax_transmit>; 

defines  the  maximum  number  of  nodes  in  the  shading 
tree,  the  minimum  distances  that  a  valid  ray  inter¬ 
section  must  be  from  the  origin  of  a  ray  to  prevent 
false  cay  bounces  and  false  shadows  caused  by 
round-off  errors,  and  linear  intensity  attenuation 
for  reflected  and  refracted  rays. 
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IMAGE_ OUTPUT:  <directory_name>,  <file_name>; 

creates  a  file  where  the  generated  raster  image  will 
be  written  by  the  EXECUTE  command. 

POINT_OUTPUT:  <di rectory_name>,  <file_name>; 

creates  a  file  where  the  image  sample  points  will  be 
written  by  the  EXECUTE  command.  This  file  can  be 
later  converted  by  a  post-processor  into  a  raster 
image.  Post-processor  FOCUS  generates  images  which 
are  focused  and  have  a  depth  of  field;  post-proc¬ 
essor  BLUR  adds  motion  blur. 

PAINT_TABLE:  <identif ier>,  l<from_x>,  <from_y>,  <from_z>, 
<from_red>,  <f rom_green>,  <from_blue>, 

<to_x>,  <tO_y>r  <to_z># 

<to_red>,  <to_green>,  <to_blue>]; 
defines  a  paint  table  which  contains  a  sequence  of 
3D  vectors  specified  by  their  initial  and  final 
points  and  corresponding  intensity  values.  The 
intensity  values  are  linearly  interpolated  along  the 
direction  of  the  vector  and  used  as  ambient 
intensity  to  be  painted  on  a  surface. 

INTENSITY_MAP:  <identif ier >,  <di rectory_name>,  <file_name>; 

defines  a  color  image  which  is  to  be  used  as  an 
ambient  intensity  map  of  a  color. 

TEXTURE_MAP:  <identi f ier >,  <di rectory_name>r  <file_name>; 

defines  a  monochrome  image  which  is  to  be  used  as  a 
texture  function  of  a  color. 


236 


COLOR:  <identif ier>,  <color_coef f icients>  &  <image_map>, 


<image_map_coef f icients>  &  <texture_map>, 
<texture_map_coef ficients>  &  <paint_table>; 
defines  the  color  properties  of  a  surface  element. 
This  command  is  identical  to  the  COLOR  command  in 
SOED  but  the  source  of  ambient  intensity  can  also  be 
a  color  image  or  a  paint  table,  and  a  texture  func¬ 
tion  defined  in  a  monochrome  image  may  also  be 
added. 

LIGHT:  <identif ier>,  <x>,  <y>,  <z>,  <red>,  <green>,  <blue>, 
{shadow  I  noshadow) ,  <contrast>; 

defines  the  name,  position,  color  and  intensity  of  a 
point-light  source  which  illuminates  the  scene. 
Optionally,  shadows  cast  by  the  light  source  can  be 
computed  with  a  specified  contrast. 

LOAD:  <directory_name>,  <file_name>; 

loads  a  structured  surface  file  generated  by  the 
SUED  processor. 

PATCH_TREES:  <identif ier >,  <color>,  <di rectory_name>, 
<file_name>; 

loads  a  file  of  structured  bicubic-patch  quadtrees 
from  a  file. 

SPHERE:  <parameter_list>; 

VERTEX:  <parameter_list>; 

PLANE:  <parameter_list>; 

POLYHEDRON:  <pararoeter_list>; 


237 


PATCH:  <parameter_list>; 

QUADRIC_COEFFS;  <parameter_list>; 

QOADRIC_BOUNDS :  <parameter__list>; 

enter  unstructured  surface  elements  as  in  the  SUED 
processor. 

STATUS:  <what>; 

lists  the  current  status  and  contents  of  the 
specified  data  in  the  processor. 

LIST:  <what>; 

defines  the  parts  of  STRAW  data  that  are  to  be 
listed  during  image  generation. 

INITIALIZE; 

initializes  the  STRAW  data  base. 

EXECUTE; 

computes  a  synthetic  image  and  writes  it  to  the 
currently  opened  output  image  and  sample  files.  In 
an  animated  sequence  any  part  of  the  data  base  can 
be  modified  between  successive  EXECUTE  commands. 

EXIT; 

terminates  execution  of  the  processor. 

STEREO:  (left  I  right) ,  <view_x>,  <view_y>,  <view_z>,  <angle>; 

replaces  the  current  camera  parameters  by  those  of  a 
complementary  stereo  view  computed  from  the 
specified  view  point  and  angle  of  separation. 

INPUT:  <file_name>; 

directs  the  command  processor  to  read  subsequent 
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commands  from  the  specified  file. 


HELP; 


lists  the  currently  implemented  commands  and  their 
syntax. 


•I 


m 


Unclassified 


«CUWTY  CLASSIFICATION  or  THIS  PAGE  (When  D«re  Enrered) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


1.  REPORT  number 

IPL-TR-033 


2.  GOVT  ACCESSION  NO.  3.  RECIPIENT'S  CATALOG  NUMBER 


tkkJkllAitL 


4.  TITLE  (e nd  Submit) 


Generating  Three-Dimensional  Surface  Models 
of  Solid  Objects  From  Multiple  Projections 


S.  TYPE  OF  REPORT  b  PERIOD  COVEREO 

Technical 


<-  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHOR^ 

Michael  Potmesil 


8.  CONTRACT  OR  GRANT  NUMBERfiJ 

N00014-82-K-0301 


8.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Image  Processing  Laboratory 
Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 


10.  program  element,  project,  task 

AREA  A  WORK  UNIT  NUMBERS 


11.  CONTROLLING  OFFICE  NAME  AND  ADORESS 

Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  VA  22217 


12.  REPORT  DATE 

October  1982 


<4.  MONITORING  AGENCY  name  a  ADDRESSfif  dll ft rent  Item  Controlling  Olllct) 


IS.  NUMBER  OF  PAGES 

239 


IS.  SECURITY  CLASS,  (ol  thlt  report; 

Unclassified 


I5a.  DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


IS.  DISTRIBUTION  STATEMENT  ( o t  thlt  Report) 


"The  United  States 'Government  is  authorized  to  reproduce  and  distribute 
reprints  for  governmental  purposes  notwithstanding  any  copyright  notices 
heron." 


17.  DISTRIBUTION  STATEMENT  (al  tht  tbtlrtel  entered  In  Block  20.  II  different  from  Rtport) 


IB.  SUPPLEMENTARY  NOTES 


IB.  KEY  WOROS  (Continue  on  reeerte  tide  It  ntcttttry  and  Identify  by  block  number •) 

Computer  Graphics  Image  Processing 

Computer-aided  design  Pattern  Recognition 

3D  modelling 
Solid  geometry 

20.  ABSTRACT  (Continue  on  reveree  tide  It  neceeeery  end  Identify  by  block  number) 

A  recurring  problem  in  computer  vision  and  related  fields  is  that  of  generating 
computer  models  of  physical  objects.  This  thesis  presents  a  method  for  con¬ 
structing  such  models  in  the  form  of  three-dimensional  surface  or  volume  descrip¬ 
tions.  The  surface  models  are  composed  of  curved,  topologically  rectangular, 
parametric  patches.  The  data  required  to  define  these  patches  are  obtained  from 
muntple  photographic  images  of  the  object  illuminated  by  a  rectangular  pattern 
of  lines.  The  projection  of  the  pattern  on  the  surface  of  the  object  (over) 


Unclassified _ 

SECURITY  CLASSIFICATION  or  THIS  PAOEfWii  Data  Zntand) 

traces  curves  which  define  the  boundaries  of  the  patches.  The  3D  description 
of  the  patches.  The  3D  description  of  the  patches  is  reconstructed  by  photo- 
grarrmetric  techniques  from  two  or  more  images  of  the  projected  pattern.  A 
calibration  stand,  in  which  the  object  is  placed,  permits  determination  of  the 
camera  geometry  directly  from  image  data. 

This  method  generated  3D  surface  descriptions  of  only  those  parts  of  the 
object  that  are  illuminated  by  the  projected  pattern,  and  also  are  visible 
in  at  least  two  images.  A  complete  model  of  the  object  is  obtained  by  re¬ 
peating  this  reconstruction  process  for  various  arbitrary  orientations  of 
the  object  until  descriptions  covering  the  entire  surface  have  been  obtained. 
Since  each  description  is  defined  in  its  own  3D  object  space  and  2D  parameter 
space,  a  surface-matching  procedure  is  developed  to  register  spatially  the 
surface  descriptions  in  a  common  object  space.  This  procedure  searches  for 
a  3D  rigid  transformation  of  the  surface  descriptions  which  minmizes  their 
shape  difference.  Once  the  surface  descriptions  are  in  the  same  object  space, 
they  are  also  merged  into  a  common  parameter  space.  This  match-and-merge 
process  is  iteratively  repeated  for  pairs  of  surface  descriptions  until  a 
completer  model  of  the  object  is  assembled. 


SECURITY  CLASSiriCATIO*  OF  TUI*  PAGErWiw  .:>•»«  Enttrtd) 


